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1  Introduction 


Inverse  scattering  has  been  an  active  field  of  research  in  science,  mathematics,  and  engi¬ 
neering  over  the  past  several  decades  (see  e.g.  [2],  [3],  [5],  [7],  [9],  [10],  [11],  [13]).  It 
has  applications  in  a  wide  range  of  fields,  such  as  radar,  medical  imaging,  oil  exploration, 
microscopy,  etc.  However,  constructing  reliable  and  efficient  algorithms  for  the  solution  of  in¬ 
verse  scattering  problems  involves  three  major  difficulties.  First,  inverse  scattering  problems 
are  highly  nonlinear  except  for  the  one-dimensional  case,  where  the  nonlinear  problem  can 
be  reduced  to  a  linear  one,  although  the  procedure  of  the  reduction  is  numerically  unstable. 
Second,  numerical  stability  tends  to  be  a  problem,  except  in  one  dimension.  Third,  the  time 
and  memory  requirements  are  beyond  the  capabilities  of  present  computers  via  currently 
used  methods.  In  this  paper,  we  construct  numerical  algorithms  for  the  solution  of  inverse 
scattering  problems  in  the  acoustic  environment  in  three  dimensions.  Our  inverse  scattering 
scheme  assumes  that  the  speed  c{x,  y,  z)  of  propagation  of  sound,  the  density  p(x,  y,  z)  and 
the  attenuation  '-y{x,y,z)  are  independent  of  the  variables  x,y,  so  that  c{x,y,z)  =  c{z), 
p{x,y,z)  =  p{z),  'y{x,y,z)  =  'y{z)]  an  acoustic  medium  possessing  these  properties  will  be 
referred  to  as  a  layered  medium,  or  layered  environment. 

1.1  Statement  of  the  Problem 

The  inverse  scattering  problem  is  the  problem  of  reconstructing  the  various  parameters 
of  scattering  objects,  such  as  the  density,  the  speed  of  sound,  and  the  attenuation,  with 
the  knowledge  of  the  incident  and  the  scattered  field.  Below  is  the  formal  mathematical 
formulation  of  the  three-dimensional  inverse  scattering  problem  in  a  layered  acoustic  medium. 

1.1.1  The  Helmholtz  Equations 

The  inverse  scattering  problem  we  investigate  arises  from  the  time  domain  wave  equation 

^2  I 

— t)  =  c^(x)  ■  p(x)  V  ■  (-^V^(x,  t)),  (1) 

where  'ip{x,t)  is  the  value  of  the  scalar  field  at  a  point  x  at  time  t,  c(x)  is  the  local  speed  of 
wave  propagation  at  a  point  x,  and  p(x)  is  the  density  at  a  point  x.  In  order  to  solve  (1), 
we  assume 

(2) 

where  k  is  a  complex  number  with  non- negative  imaginary  part,  and  Cq  is  the  speed  of  wave 
propagation  outside  of  the  scattering  structure.  Substituting  (2)  into  (1),  we  obtain 

P(^)  ^  +  k^  ■  =  0.  (3) 

P[X)  C^(X) 

Since  the  inverse  scattering  scheme  assumes  that  the  speed  c(x,  y,  z)  of  propagation  of  sound, 
the  density  p{x,  y,  z)  and  the  attenuation  7(0;,  y,  z)  are  independent  of  the  variables  x,  y,  i.e., 
c{x,y,z)  =  c{z),  p{x,y,z)  =  p{z),  'y(x,y,z)  =  7(2:),  (3)  can  be  rewritten  by  the  formula 

_  ^  ^  2  .  eg 
)  dz  dz  c^{z 


y  •  =  0.  (4) 


(9Vfc _ i_ 

dx"^  dz"^  p{z 
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Throughout  this  paper,  we  use  the  notation 


^— =  1 +g(^) +i-7(^),  (5) 

c^{z) 

where  g(z)  and  j(z)  are  known  as  potential  and  attenuation  of  the  layered  acoustic  medium, 
and  that  p,  q,  ■y  E  Cg([0, 1]),  i.e.,  p,  q,  7  are  twice  continuously  differentiable  everywhere,  and 
are  dehned  by  the  formulae 


p{x)  =  p(0)  =  pi,  for  all  a;  <  0,  (6) 

p(a;)  =  p(l)  =  p2,  for  all  a;  >  1,  (7) 

q(x)  =  q(0)  =  gi,  for  all  a;  <  0,  (8) 

q(x)  =  g(l)  =  q2,  for  all  a;  >  1,  (9) 

y'(x)  =  7(0)  =  7i,  for  all  a;  <  0,  (10) 

7(a;)  =  7(1)  =  72,  for  all  a;  >  1.  (11) 


Suppose  now  that  the  angle  of  incidence  with  respect  to  the  normal  to  the  x-y  plane  is 
6*,  and 


^k{x,z) 

Substituting  (12)  into  (4),  we  obtain 


=  e 


ikx  sinO 


(l>{z). 


(12) 


(j)''{x,  k) 


Ax) 

p{x) 


■  (j)'{x,  k)  +  k'^  ■  {1  +  q{x)  +  i  ■  y'(x) 


A)  ■  (j){x,  k)  =  0, 


(13) 


where 

a  =  sin{6).  (14) 

Equation  (13)  is  the  well-known  scalar  Helmholtz  equation  in  layered  acoustic  media.  For 
any  complex  k,  we  consider  solutions  of  the  Helmholtz  equation  0+(a;,  k)  and  0_(a;,  k)  dehned 
by  the  formulae 

0+(a;,  k)  =  (l)inc+{x,  k)  +  (pscat+ix,  k),  (15) 

(f)-{x,  k)  =  k)  +  (f)scat-{x,  k) ,  (16) 

with 

(j)inc+iX:  k)  =  for  all  a;  <  1,  (17) 

(j)inc+{x,  k)  =  72-^2  ^ 

(j)inc-{x,  k)  =  for  fol  x  <  0,  (19) 

^inc-{x,  k)  =  e-*^Vl+'?2+i72-a2x^  ^  ^  ^20) 

and  0scat+,  4>scat-  satisfying  the  outgoing  radiation  boundary  conditions 

Aat±i^:  k)+ik  v^l  +  gi  +*71  (j)scat±i0,  k)  =  0,  (21) 
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(22) 


Combining  equations  (13)  -  (20),  we  obtain  the  equations 

<P'scat+{x,  k)  -  •  (I)'scat+{X,  /c)  +  ■  (1  +  q{x)  +  i  •  -^{x)  -  O^)  ■  (t)scat+{x,  k) 

=  -{k^  {{q  -qi)+i{l-  7i))  -^ik  +  * 7i  -  «')  •  e*' 

p[x) 


for  all  a;  <  1,  (23) 

4>"scat+{x,  k)  -  /c)  +  ■  (1  +  g(a;)  +  i  ■  7(0;)  -  a^)  ■  (pscat+ix,  k) 

=  -{k‘^{{q-q2)  +*(7-72))  -  a/1  +  ^2  +  *72  -  a^)  •  e*fcVi+52+» 72-^2 x 

p(a;) 

for  all  a;  >  1,  (24) 

(x) 

0"cat-(a^,  ■  <Pscat-ix,  k)  +  k‘^  ■{!  +  q{x)  +  i  ■  j{x)  -  a^)  ■  (pscat-i^,  k) 

_  flJ2  ^  I  „•  I  P  (^)  „•  L  /i  i  I  I  2.,  l+Oi  +i  71  -oi^  x 


=  -{k^  {{q  -  gi)  +  *  (7  -  7i))  +  ^^k  v/l  +  gi+*7i-«')  ■  g- 

p(a;) 

for  all  a;  <  0,  (25) 


<P'Lat-ix,  k)  -  •  (l)'scat-{x,  k)  +  k‘^  ■  {I  +  q{x)  +  i  ■  7(a;)  -  a^)  ■  (t)scat-{x,  k) 

p[x) 

=  -{k^  {{q  -  g2)  +  p7  -  72))  +  i  k 

p[x) 

for  all  a;  >  0.  (26) 

Combining  the  equations  (23)  -  (26)  with  the  equations  (6)  ~  (11),  we  observe  that,  for  any 
complex  k,  there  exist  complex  numbers  pi±{k)  and  po±{k)  such  that 

4>scat±{x,  k)  =  p^±{k)  ■  ^11  ^  (27) 

(t>scat±{.x,  k)  =  po±{k)  ■  for  all  a;  <  0;  (28) 

combining  (18),  (19),  (27),  and  (28),  we  obtain 

0+(a;,  k)  =  {l+  pi+{k))  ■  Vi+^^+i 72-^2 ^  ^  (29) 

(t)-{x,k)  =  (l+/ao-(A;))-g-'^v^^+^^+^^^^foralla;<0.  (30) 

Thus,  for  any  complex  k,  the  boundary  value  problems  for  0+,  0_  (equations  (13)  -(22))  are 
reformulated  as  initial  value  problems  (equations  (13),  (29),  and  (30)).  Furthermore,  for  any 
k  E  C,  coefficients  1  -|-  pi+{k)  and  1  -|-  po-{k)  in  (29),  (30)  are  both  nonzero.  For  example, 
if  1  -|-  pi+{k)  =  0,  then  0+(l,  k)  =  0+(l,  k)  =  0,  thus,  due  to  uniqueness  theorem  on  initial 
value  problems,  0+(a;,  /c)  =  0  for  all  x  E  R,  i.e., 

(Pscat+ix,  k)  =  -<t>,nc+{x,  k)  =  _ g*  ^  x ^  (3I) 

for  all  a;  <  1,  which  contradicts  (28).  Similarly,  we  can  prove  that  l+po-{k)  is  also  non- zero. 
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1.1.2  The  Impedance  Functions 

Denote  the  upper  half  complex  plane  by  .  For  any  k  G  we  dehne  the  impedance 
functions  p+,P-  :  (-R,  C^)  C  hy  the  formulae 


p+[x,k)  =  .  ^  .  ,  V 

i  k  p[x)  (p+yX-,  k) 

(32) 

/  4>'-{x,k) 

P-{x,k)  =  /  N  ,  /  ,  V 

—i  k p[x)  (p-{x,  k) 

(33) 

where  0+,  are  solutions  of  the  Helmholtz  equation  (13)  dehned  by  (18)  and  (19),  p  is  the 
density  of  the  scattering  structure. 

Remark  1.1  It  is  easy  to  see  that  the  impedance  functions  p+,  p_  are  independent  of  the 
nonzero  coefficients  l  +  /ii+(A;)  and  l+p,o-{k)  in  initial  conditions  (29),  (30).  Therefore,  for 
simplicity,  the  initial  conditions  (29),  (30)  are  reformulated  as 

k)  =  ^  (34^ 

k)  =  for  all  a;  <  0.  (35) 

The  solutions  0+,  as  solutions  of  equation  (13)  subject  to  boundary  conditions  (34),  (35) 
only  differ  from  those  subject  to  boundary  conditions  (29),  (30)  by  constants. 

Therefore,  the  inverse  scattering  problem  for  the  equation  (13)  is  stated  as  follows: 

Suppose  that  the  impedance  function  p+(0,  k)  is  known  for  appropriately  chosen  frequen¬ 
cies  {kj,j  =  1,  2, ...,  M},  reconstruct  the  density  p{z),  the  potential  q{z),  and  the  attenuation 
y{z),  in  the  interval  [0, 1]  with  the  error  that  rapidly  decreases  with  increasing  M . 

This  paper  is  devoted  to  the  construction  of  an  algorithm  for  the  solutions  of  the  above 
problem. 

1.2  Overview 

As  discussed  in  [3],  four  types  of  algorithms  exist  for  the  solution  of  inverse  scattering 
problems. 

1.  Linearized  inversion  schemes,  which  approximate  the  inverse  scattering  problem  by 
the  problem  of  inverting  the  appropriately  chosen  linear  operator  (see  e.g.  [9]). 

2.  Methods  based  on  nonlinear  optimization  techniques,  which  obtain  the  scattering 
parameters  iteratively  by  solving  a  sequence  of  forward  scattering  problems  (see  e.g.  [13]). 

3.  Gel’fand-Levitan  and  related  techniques,  converting  the  Helmholtz  equation  into  the 
Schrodinger  equation,  and  solving  an  inverse  problem  in  the  form  of  a  linear  Volterra  integral 
equation  for  the  latter  (see  e.g.  [7]). 

4.  Methods  based  on  trace  formulae,  which  connect  the  high-frequency  behavior  of 
the  solutions  of  the  Helmholtz  equation  to  the  parameters  of  the  scattering  objects  to  be 
recovered  (see  e.g.  [11]). 

In  this  paper,  we  introduce  an  algorithm  for  the  solution  of  inverse  scattering  problems 
in  layered  acoustic  media.  The  procedure  falls  into  the  category  4,  and  is  an  extension  of 
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the  procedure  of  [3],  in  the  sense  that  while  the  latter  recovers  the  parameters  of  a  layered 
medium  in  which  only  the  speed  of  sonnd  is  permitted  to  change,  the  algorithm  of  my  thesis 
assumes  that  the  speed  of  sound,  the  density,  and  the  attenuation  are  variable,  and  recovers 
all  three. 

The  inverse  scattering  schemes  we  constrnct  are  based  on  a  collection  of  so-called  trace 
formulae,  and  can  be  viewed  as  extension  of  the  work  started  in  [3],  where  the  observation  is 
made  that  (at  least  in  layered  media)  it  is  possible  to  constrnct  inverse  scattering  algorithms 
that,  given  a  smoothly  varying  medinm,  reqnire  few  measnrements  to  reconstruct  it.  More 
specihcally,  given  a  medium  whose  parameters  c,  p,  and  7  have  m  >  1  continuous  derivatives, 
and  data  measured  for  all  frequencies  u  on  the  interval  [—a,  a],  the  error  of  the  reconstrnction 
decays  as  as  a  — cx).  In  this  respect,  the  algorithm  of  [3]  is  similar  to  the  Fourier 

Transform,  and  a  strong  argnment  is  made  that  this  is  a  very  desirable  property.  While 
the  algorithm  of  [3]  assumes  that  the  parameters  p  and  7  are  constant  and  the  parameter 
c  depends  on  2:,  the  schemes  of  this  paper  reconstrnct  c,  p,  and  7,  provided  that  they  only 
depend  on  the  coordinate  2:. 

The  paper  is  organized  as  follows.  In  Section  2,  we  snmmarize  several  well-known  math¬ 
ematical  facts  to  be  used  in  the  paper.  In  Section  3,  we  introduce  analytical  tools  to  be 
used  in  the  construction  of  the  algorithm.  Section  4  states  the  algorithm  in  details,  and 
a  complexity  analysis  is  included.  In  Section  5,  several  numerical  examples  are  used  to  il¬ 
lustrate  the  performance  of  the  algorithm.  Finally,  Section  6  contains  generalizations  and 
conclusions. 


2  Analytical  Preliminaries 

In  this  section,  we  summarize  several  well-known  mathematical  facts  to  be  used  in  the 
sections  below.  These  facts  are  given  withont  proofs,  since  all  of  these  are  either  found 
in  [1],  [4],  and  [12],  or  easily  derived  from  well-known  results. 

2.1  Notation 

In  this  paper,  we  denote  the  upper  half  complex  plane  by  . 

For  any  function  f  :  R  ^  R,  f  E  c™([0, 1])  states  that  is  continnous  everywhere, 
and  that  f{x)  =  fi  for  all  a;  <  0,  f{x)  =  f2  for  all  a;  >  1,  where  fi  =  /(O),  /2  =  /(I). 
Snppose  that  for  any  a  >  0,  the  region  K{a)  C  C  is  dehned  by  the  formula 

K{a)  =  {k\k  G  C,Im{k)  >  0,  \k\  >  a}.  (36) 

2.2  Basic  Lemmas 

In  this  section,  we  introdnce  several  basic  lemmas  to  be  used  in  the  sections  below  following 
closely  to  [3]. 

For  a  fnnction  /  :  [a,  b]  R^  and  integer  n  >  2,  the  n-point  trapezoidal  rule  Tn  is  dehned 
by  the  formula 

+  (37) 
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where 


(38) 


h  = 


b  —  a 


n  —  1  ’ 

and  is  widely  used  as  an  approximation  to  the  integral  J^f{x)dx.  It  is  second  order  con¬ 
vergent,  as  long  as  the  function  /  has  two  continuous  derivatives  on  [a,  b].  In  other  words,  if 
/  G  C^[a,  6],  then 


f(x)  dx  =  r„(/)  +  0(A=). 


(39) 


Lemma  2.1  Suppose  that  f  G  c\R)  with  I  a  nonnegative  integer.  Suppose  further  that 
=  0  for  0  <  j  <  I,  f^^'>  is  absolutely  eontinuous.  Then  there  exists  a  positive  number 
c  such  that 

fX  1  j  1  «  +  i 

/  )  ;,-i(j,)  +  (  )  (40) 


J  =  1 


with  b  :  R  X  C  an  absolutely  continuous  function  of  x  E  [0, 1]  such  that 

\b{x,k)\<c,  (41) 

for  all  X  G  [0, 1],  k  G  .  Furthermore,  if  f{x)  =  0  for  all  x  >  D  with  D  a  positive  number, 
then 

\b{x,k)\<c,  (42) 

for  all  {x,  k)  E  R  X  . 

Lemma  2.2  Suppose  that  a  :  [0, 1]  — i?  and  b  :  [0, 1]  — C  are  two  continuous  functions, 

and  that  a(x)  >  0,  for  all  x  G  [0, 1].  Then  for  any  two  solutions  u  and  v  of  the  second  order 

ODE 

{a{x)  (j)'{x)y  +  b{x)  (f^x)  =  0,  (43) 

there  exists  a  constant  c  such  that 


a{x){u{x)v'  (yx)  —  v{x)u{x))  =  c 


(44) 


for  all  X  E  [0, 1].  Furthermore,  c  0  if  and  only  if  u  and  v  are  linearly  independent.  (The 
expression  W{u,v)  =  u{x)v\x)  —  v\x)u{x)  is  referred  to  as  the  Wronskian  of  the  pair  u  , 
v). 

Lemma  2.3  (GronwalTs  ineguality)  Suppose  that  u,v,w  :  [0,  a]  — R  are  three  continuous 
and  nonnegative  functions,  satisfying  the  ineguality 


w{x)  <  u{x)  -1-  /  v{t)  w(t)  dt 


(45) 


for  all  X  E  [0,  a].  Then, 


w{x)  <  u{x)  +  /  u{t)v{t)e-^t  dt, 


(46) 


for  all  a;  G  [0,  a  . 


Lemma  2.4  Suppose  that  a  :  C  ^  C  is  an  entire  funetion  and  that  A  :  R  x  C  ^  (jnxn 
an  n  X  n-matrix  whose  entries  aij{x,z),  i,j=l,...,n  are  continuous  funetions  of  x  and  entire 
functions  of  z  for  all  x  &  R  Then  for  any  z  E  C,  the  differential  equation 

y\x,z)  =  A{x,z)-y{x,z)  (47) 


subject  to  the  initial  condition 


y(0)  =  c(z) 


(48) 


has  an  unique  solution  y{x,  z)  for  all  x  E  R.  Moreover,  y{x,  z)  is  an  entire  function  of  z. 


2.3  Schrodinger  Equation  and  Riccati  Equation 


This  section  describes  the  basic  facts  about  the  Helmholtz  equation  and  its  connections  with 
the  Schrodinger  equation  and  the  Riccati  equation  in  the  context  of  scattering  problem. 
Lemma  2.5  describes  the  fact  that  a  Schrodinger  equation  with  outgoing  radiation  conditions 
can  be  converted  into  a  second  kind  integral  equation  with  the  Green’s  function  of  the 
corresponding  Helmholtz  equation.  Lemma  2.6  gives  the  mathematical  form  of  the  Green’s 
function  for  Helmholtz  equation  with  outgoing  radiation  conditions. 


Lemma  2.5  Suppose  that  Gk  ■  [0, 1]  x  [0, 1]  — >  G  is  the  Green’s  function  of  the  boundary 
value  problem 

^jJ"{x,k)  +  k'^  ^jJ{x,k)  =  0  (49) 

f'{0,k)+ikij{0,k)  =  0  (50) 

^jJ'{l,k)  —  ik^jJ{l,k)  =  0  (51) 

for  any  complex  k  ^  0.  Then  the  boundary  value  problem 


ip''{x,  k)  +  (A;^  +  rj{x))  ip{x,  k)  =  f{x,  k) 

k)  +  i  /c'^(0,  k)  =  0 
iffl,  k)  —  ikipfi,  k)  =  t) 
is  equivalent  to  a  second  kind  integral  equation 

ip{x,  k)  =  —  Gk{x,  f)  r]{f)  'ijj{t,  k)  dt  +  g{x,  k) 

Jo 

with  f,g:[0,l]xG^G  and  g  defined  by  the  formula 

g{x,k)=f  Gk{x,t)  f{t,k)dt. 

Jo 

Lemma  2.6  For  any  complex  k  ^  0,  the  Helmholtz  equation 


(52) 

(53) 

(54) 

(55) 


(56) 


fi"{x,  k)  +  'ip{x,  k)  =  t) 


(57) 


with  the  outgoing  radiation  conditions  (21),  (22)  has  the  Green’s  function 


Gk{x,t) 


1  (  Q^k{t  x)  X  <t, 

2ik\  X  >t. 


(58) 
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The  following  lemma  connects  the  solutions  of  the  Helmholtz  equation  to  those  of  Schrodinger 
equation  via  direct  transform. 

Lemma  2.7  Suppose  that  qj'y,  p  :  R  ^  R  are  c^-functions  sueh  that  1  +  q{x)  —  >  0, 

7(0;)  >  0  ,  p{x)  >  0,  for  all  x  E  R.  Suppose  further  that  the  functions  n,t  :  R  ^  C , 
p,  g  :  C  ^  C  are  defined  by  the  formulae 


n{x)  =  a/I  -I-  q{x)  +i7{x)  —  a'^, 

(59) 

px 

t{x)  =  /  n(r)  dr, 

Jo 

(60) 

V{i)  =  7(1  +  ^  ( (2  -  3  ■  (1  +  q(x)  +  i  j(x)  -  a^f 

4  V  p(x)  p(x) 

-(q"(x)  +ij"(x))  ■  (1  +  q(x)  +ip(x)  -  a^)  +  ^(q(x)  +i7'(x)f^  (61) 

g(t)  =  f(x)  ■  p~^(x)  ■  (1  +  q(x)  +  i7(x)  —  a^)  .  (62) 

Finally,  suppose  that  the  function  (f)  :  R  x  C  ^  C  satisfies  the  equation 

4>''{x,  k)  —  ^  (j)'{x,  k)  +  (1  +  q{x)  +  f  7  —  a^)  (j){x,  k)  =  f{x)  (63) 

p{x) 

and  the  function  ip  :  C  x  C  ^  C  is  defined  by  the  formula 

fj{t,  k)  =  p~^{x)  ■  (1  +  q{x)  +  i7{x)  —  ■  (j){x,  k).  (64) 

Then  the  function  fj  satisfies  the  Schrodinger  equation 

%l)''(t,  k)  +  {k^  +  p{t))  •  ij{t,  k)  =  g{t).  (65) 

Remark  2.1  Lemma  2.7  provides  a  connection  between  the  solutions  of  the  Helmholtz  equa¬ 
tion  (63)  and  those  of  the  appropriately  chosen  Schrodinger  equation  (65).  This  connection 
will  be  used  in  the  following  chapter  as  an  analytical  tool.  However,  it  is  not  useful  in  nu¬ 
merical  computations  since  the  connection  between  p  and  q  (  see  equation  (61))  is  generally 
ill-conditioned. 

Corollary  2.8  Suppose  that  under  the  conditions  of  the  preceding  lemma  that  q,7,p  G 
Co([0, 1]).  Suppose  further  that  the  functions  ip-  ■.  C  xC  ^  C  are  defined  by  the  formulae 

i/j+it,  k)  =  p~^{x)  ■  (1  -h  q{x)  +i7{x)  -  a^Y  ■  0+(a;,  k),  (66) 

ip-{t,  k)  =  p~^{x)  ■  (1  -I-  q{x)  +  ipfx)  —  k).  (67) 

Then  satisfy  the  ODEs 

^+(t,  k)  +  (/c^  +  p{t))  •  fj+Y,  k)  =  0,  (68) 
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(69) 


k)  +  (/c^  +  'r]{t))  ■  'ip-it,  k) 


subject  to  the  boundary  conditions 


0, 


for  all  Re(t)  >  Re(Ti),  and 


^/J-{t,  k)  =  p^^l  +  qi-a  +  i7i)4  e 


T  k t 


for  all  Re(t)  <  0  with  Ti,  7^  0  defined  by  the  formulae 

Ti=  a/1  +  q{x)  +  i  'y{x)  —  a‘^  dx 

Jo 

^{k)  =  p“^  (1  +  g2  -  a'  +  *72)3 

Furthermore, 

p+{x,k)  =  Vl  +  q{x)  +  i  7(3;)  -a^  .  + 

i  k  p[x)  ip+\t,  k) 

i  ®  “  i  (1  +  ■  (^'(^)  +  * 

i  k  p{x) 


(70) 

(71) 


(72) 

(73) 


(74) 


P-{x,  k)  =  a/1  +  q{x)  +  i 7(0;)  — 

-  i(l  +  ^(^)  -  tt 


f)'_{t,k) 


i  k  p{x)  if-{t,  k) 


i  k  p{x) 


(75) 


Observation  2.2  Suppose  that  q,'-^,  p  &  Cg([0, 1]).  Then  according  to  Lemma  2.1  and  Corol¬ 
lary  2. 8, 

t  =  a/1  +  gi  +  i7i  -  a^a;  (76) 

and  consequently 

(t)+{x,  k)  =  ■  (1  +  (?i  +  i  7i  -  q;^)“5  •  ^7+(t,  k),  (77) 


/or  allx  <  0.  Now,  suppose  the  function  is  defined  by  formula  (66).  Defining  the  scattered 
field  'ipscat+  '.  C  X  C  ^  C  by  the  formula 

fj+{t,  k)  =  p^^  ■  (a/1  +  gi  +i7i  -  a2)2  ■  (e*^‘  +  fjscat+it,  k)),  (78) 


we  immediately  obtain  the  Schrodinger  equation 


J^'Lt+{t^k)  +  {k‘^ +  p{t))f)scat+{t,k)=  p{x)  2  .  (1 +  g(a;) +i7(a;)  -  a^)  "■ 

((g  -  gi)  +  ^  (7  -  7O)  +  ^lk  v/l  +  gi+^7i-a2)  .  ^  ^79^ 

p{x) 

subject  to  outgoing  radiation  conditions  (21),  (22). 
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The  following  lemma  introduces  the  Riccati  equations  satisfied  by  the  impedance  functions 
P+,  P-.  They  are  obtained  from  the  dehnitions  of  impedance  functions  (32)  (33)  and  the 
Helmholtz  equation  (13). 

Lemma  2.9  Suppose  that  under  the  eonditions  of  the  preceding  lemma, 


'f+{xo,ko)  7^  0, 

(80) 

fj-{xo,ko)  0, 

(81) 

at  some  point  {xq,  ko)  E  Rx  C.  Then  there  exists  a  neighborhood  D  of  (xq,  ko) 
impedance  functions  p+,p-  satisfy  the  Riccati  equations 

such  that  the 

/  r  1\  -1  r2r  1\  1  +  g(a^)  +  *  7(a^)  - 

p+{x,k)=  ikp{x)-{p^[x,k)  ), 

p^[x) 

(82) 

!  t  1\  •  /  /  \  /  2  /  /  N  1  +  +  i  l{x)  -  , 

p_{x,k)=ikp{x)-{pf{x,k) 

P  \x) 

(83) 

for  all  {x,  k)  E  D. 

Observation  2.3  Combining  formulae  (34),  (35),  we  easily  observe  that 

p+(x,«:)=  +  forallx>l, 

P2 

(84) 

P.{x,k)=  f„rallx<0, 

(85) 

for  all  complex  k  ^  0. 

Observation  2.4  If  'y{x)  =  0  for  all  x  &  R,  it  is  easy  to  see  from  equations  (82),  (83), 
(84),  and  (85)  that 

p+{x,k)  =p+{x,-k),  (86) 

p^{x,  k)  =  p^{x, —k).  (87) 

for  all  real  k,  since  p+{x,  k)  and  p+{x,  —k)  satisfy  identical  differential  equations  and  bound¬ 
ary  conditions,  and  the  same  is  true  for  p-{x,  k)  and  p-(x,  —k),  too. 


3  Mathematical  Apparatus 


In  this  section,  we  introduce  analytical  tools  to  be  used  in  the  construction  of  the  algorithms 
of  this  paper.  This  section  discusses  the  following  three  facts. 

(A)  For  any  x  E  R,  the  impedance  functions  p+{x,  k),  p-{x,  k),  dehned  by  (32),  (33),  are 
analytic  functions  of  k  in  the  upper  half  plane  .  Furthermore, 


V+{.x,k)  =  — - 

p{x) 

1 

i  k 


a/1  +  q{x)  +  iy  -  q;2 

p{x)  ■  {q\x)  +  i  Y(x))  -  2  ■  (1  +  g(x)  +  i  ^(x)  -  a^)  ■  p'(x) 

4-p^{x)-{l  +  q{x)+i'y{x)-a^)  ^  ^  ’ 


12 


k) 


^  ■  a/1  +  q{x)  +i7  -  a2 
[x) 

]_  p{x)  ■  {q\x)  +  i  i{x))  -  2  ■  (1  +  q{x)  +  i  7(3;)  -  a^)  ■  /(x)  p.,  2a  .oqA 

4-p2(a;).(i  +  g(a;)+,q,(a;)_a2)  ^  ^  I  J 

as  |/c|  — cx)  for  all  a;  G  -R,  /c  G  C"*". 

(B)  For  any  fixed  x  E  R,  the  difference  between  the  impedance  functions  p+{x,  —k)  and 
P-{x,k)  decays  like  a  constant  times  k~"^,  where  k  E  R,  and  m  is  the  smoothness  of  the 
scatterer.  In  other  words, 

p+{x, —k)  —  p_{x,k)  =  0{k~'^),  (90) 

as  \k\  00,  k  E  R. 

(C)  For  any  a  >  0,  and  all  x  E  R, 

p{x)  ■  {q'{x)  +  i  7/a;))  —  2  ■  p'{x)  ■  (1  +  q{x)  +  i  ^(x)  —  a^) 

=  — {1  +  q{x)  +  i'^{x)  —  a^)  p^{x)  [  {p+{x,k)  —  p-{x,k))  dk  +  (91) 

^  J —a 

where  m  is  the  smoothness  of  the  scatterer,  p+{x,  k)  and  P-{x,  k)  are  impedance  functions 

defined  by  (32)  and  (33),  p  is  the  density  of  the  scattering  object,  q  is  the  potential,  and  7 

is  the  attenuation.  (91)  is  an  example  of  a  trace  formula. 

The  proofs  in  this  section  are  modeled  after  those  in  [3]. 

3.1  Boundedness 

This  section  establishes  the  basic  properties  of  the  impedance  functions  p_,  defined 
by  (32),  (33).  Lemma  3.1  describes  the  behavior  of  0_  in  the  vicinity  of  /c  =  0  in  the 
complex  plane.  Lemma  3.2  describes  the  properties  of  the  impedance  functions  p+,  p-  near 
k  =  0.  Lemma  3.3  shows  that  0_,  0+  are  nonzero  for  all  real  x  and  complex  A;  /  0. 

The  following  lemma  describes  the  behavior  of  0_  in  the  vicinity  of  A;  =  0  in  the 
complex  plane. 

Lemma  3.1  Suppose  that  p,q,'-^  E  Cq([0,  1]),  and  A  >  0  is  a  real  number.  Then,  there  exist 
positive  numbers  5,  e,  and  p  such  that 


\(j)+{x,k)  -  1|  <  e|A;|, 

(92) 

\(p-{x,k)  -  1|  <  £|A;|, 

(93) 

\4>+{x,  A;)  -  i  A;  a/1  +  g2  +  *72  - 

P2 

lA 

to 

(94) 

\(p'-{x,  A;)  +  i  A;  a/1  +  gi  +  i7i  -  q;2  ■ 

Pi 

iM 

VI 

(95) 

(j)+{x,k)  /  0, 

(96) 

fi-{x,k)  /  0, 

(97) 

for  all  realx  E  [—A,  A]  and  complex  k  such  that  \k\  <  5.  In  (92)-(97),  qi,q2,'Ii,'l2,  Ojre  defined 
in  Section  2.1,  a  is  define  by  (If),  (ind  0±(a;,  k)  is  the  field  at  x. 
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Proof.  Since  the  proofs  of  this  lemma  for  and  for  <p-j<p'_  are  identical,  we  only  prove 

it  for  <p_,<p'_.  Defining  two  auxiliary  functions  (pi, ip:  R  x  C~^  ^  C  hj  the  formulae 


ip{x,  k)  = 


(pi{x,  k)  =  (p-{x,  k)  -  1, 

(p'-{.x,k)  ik  /— - — - ^ 

- yv; - ^ - a/I  +  +  *7l  -  « 

p{x)  Pi 


(98) 

(99) 


and  combining  (98),  (99)  with  (13)  and  initial  conditions  (34),  (35),  we  obtain  the  linear 
first  order  ODEs 


0'i(a;,  k)  =  p(x)  ip(x,  k)  -  i  /c  a/ 1  +  gi  +  iyi  -  ^ 

Pi 

ip\x,  k)  = - --(1  +  q{x)  +i'^{x)  -  a^){(pi{x,  k)  +  1), 

p{x) 


(100) 

(101) 


subject  to  the  initial  conditions 

0i(O,A;)=O,  (102) 

ip{Q,k)  =  Q.  (103) 

We  begin  by  showing  that  there  exist  two  continuous  functions  M,  N  :  x  R^  R'^  such 

that  for  any  s  G  R'^,  M{s,t),  N{s,t)  are  monotonically  increasing  functions  of  t  and 

\(pi{x,k)\  <  M{A,\k\)\k\,  (104) 

\ip{x,k)\  <  N{A,\k\)\k\\  (105) 

Integrating  (100)  from  0  to  x,  we  have 

(pi{x,k)  =  /  (-i  /c  a/1  +  gi  +  i  7i  -  q;2  p{t) - V  p{t)ip{t,k))  dt,  (106) 

Jo  Pi 

an  substituting  into  (101)  and  integrating  the  result,  we  have 
R"  1  +  q{t)  +  i7(t)  - 


ip{x,  k)  =  —k^ 


pif) 


1+  (■ 


-i  k  a/1  +  gi  +  i  7i  -  p(r) 


Pi 


a{x,  k)  <  Cl  \kp  (  |a;|  +  -  C2  x^  \k\  +  /  {x  —  t)  a{t,  k)  dt]  , 

^  J  C\ 


+  p{T)ip{T,k))  dr  dt.  (107) 


Denoting  \ip{x,  k)\  by  a{x,  k),  we  have 

1 


where 


Cl  =  sup 

-A<x<A 


1  +  q{x)  +  i  7(0;)  — 


C2  =  sup 

-A<x<A 


p{x) 

a/1  +  gi  +  i7i  - 

Pi 


(108) 

(109) 

(110) 
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(Ill) 


C3  =  sup  \p{x)\. 

-A<x<A 

Using  Lemma  2.3,  we  rewrite  (108)  in  the  form 

a{x,k)  <\k\‘^u{x)  +  \k\^  rn(t)n(t)(l  +  e-^'^i‘^^'=' ("-*)')  ht,  (112) 

Jo 

where 

n(t)  =  Cl  |t|  +  -  Cl  C2t^  |/i;|,  (113) 

i;(t)  =  CiC3(a;-t).  (114) 

Now,  (105)  follows  immediately  with  N{A,  \  k\)  dehned  by  the  formula 

pX 

N{A,\k\)=  sup  u{x)  +  \k\^  M(t)u(t)(l  +  e-^"i"2^'(^-*)')dt  (115) 

-A<x<A  Jo 

It  is  easily  observed  from  (106)  that 

\4>i{x,  k)\  <  C2  \x\  \k\  +  c^N{A,  \k\)  \x\  \k\‘^,  (116) 

with  C2,  C3,  N{A,\k\)  dehned  by  (110),  (111),  (115),  respectively.  Therefore,  (104)  follows 
immediately  with  M{A,  \  k\)  dehned  by  the  formula 

M{A,\k\)  =  A{c2  +  c,N{A,\k\)\k\).  (117) 

Since  M{A,  t)  is  a  continuous,  monotonically  increasing  function  of  f,  there  exists  a  real  <5 
such  that 

M{A,5)-5<1  (118) 

Denoting  M{A,  5)  by  e,  N{A,  5)  by  and  observing  that  M{A,  \k\),N{A,  \k\)  are  mono¬ 
tonically  increasing  functions,  we  have 

\Mx,k)\  <M{A,\k\)\k\  <M{A,S)\k\  =e\k\,  (119) 

\i,{x,  k)\  <  N(A,  1*1)  1*1  <  iV(A  |J;|  =  A \k\\  (120) 

C3 

from  which  (93),  (95)  follow  immediately. 

Finally,  combining  (119)  with  (118),  we  obtain  (97).  □ 

The  following  lemma  describes  the  properties  of  the  impedance  functions  p_  near 
A;  =  0. 

Lemma  3.2  Suppose  that  p,q,'y  E  Co([0, 1])  and  A  >  0  is  a  real  number.  Then  there  exists 
a  real  number  <5  >  0  such  that  the  impedance  functions  p+,P-,  defined  by  (32),  (33),  are 
continuous  functions  of  {x,  k)  for  all  real  {x,  k)  G  D,  where 

D  =  {(x,  k)\x  E  [—A,  A],  k  E  C,k  ^  0,\k\  <  5}.  (121) 
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Furthermore. 


lim p+{x,  k) 

k^O 

lim  p-(x,  k) 

k^O 


a/1  +  g2  +  *  72  - 

P2 

a/1  +  gi  +i7i  - 

Pi 


(122) 

(123) 


where  qi,  q2,  7i,  72,  Pi,  P2  are  defined  in  Section  2.1,  a,  p-  are  defined  by  (I4),  (32), 
(33),  perspectively. 


Proof.  Due  to  Lemma  3.1,  there  exists  a  real  number  <5  >  0  such  that  (j)^{x,k)  7^  0, 
(f)-{x,k)  7^  0  for  all  real  {x,k)  G  D.  Therefore,  the  impedance  functions  p+,  p_  are  well- 
dehned  in  D,  and  their  continuity  follows  from  the  continuity  of  0+  ,  ,  p,  as 

well  as  their  dehnitions  (32),  (33).  (122),  (123)  are  obtained  via  the  direct  application  of 
(92)-(95)  and  (32),  (33).  □ 


Remark  3.1  Due  to  Lemma  3.2,  if  we  dehne  p+{x,  0)  =  , 

P-{x,  0)  =  then  p+,  p_  are  continuous  functions  even  at  k  =  0. 

The  following  lemma  states  that  0+,  0-,  are  nonzero  for  all  real  x  and  complex  k  ^  0. 


Lemma  3.3  For  all  x  E  R  and  complex  k  ^  0, 

<P+{x,k)fiO,  (124) 

<P'+{x,k)^0,  (125) 

<P-{x,k)fiO,  (126) 

(j)'_{x,k)  ^0.  (127) 


Proof.  Since  the  proofs  of  this  lemma  for  0+,  and  for  (p_,  are  identical,  we  only 
prove  (126)  and  (127).  We  decompose  fi-  into  two  parts  by  the  formulae 


(f)-{x,k)  =  u{x,k)  +  iv{x,k),  (128) 

4>'_{x,k)  =  u\x,k)  +  iF{x,k),  (129) 

such  that  functions  u,v:  R  x  C  ^  C  satisfy  equation  (13).  Combining  the  initial  condi¬ 
tion  (35)  with  (128),  we  obtain  that 

u{x,  k)  =  cos{k  +  qi  +  i'yi  —  x),  (130) 

v{x,  k)  =  sin(/c  \/ 1  +  qi  +  i  x) ,  (131) 

for  all  a;  <  0  and  /c  7^  0.  Lemma  2.2  states  that  the  Wronskian  W{u,v)  of  the  pair  u,  v  is 
given  by 

W{u,v)  =  a/1  -h  gi  +  i  7i  -  ■  /c,  (132) 

for  any  x  E  R.  Therefore,  for  any  A:  7^  0,  u{x,  k),  v{x,  k)  can  not  both  be  zero  simultaneously, 
nor  can  u'{x,k),  v'{x,k).  Now,  (126),  (127)  follows  immediately  given  (128),  (129).  □ 
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3.2  Asymptotics  and  Smoothness 

The  principal  purpose  of  this  section  is  to  formulate  and  prove  the  facts  (A)  and  (B)  de¬ 
scribed  in  the  beginning  of  Section  3.  We  begin  by  deriving  equations  (88),  (89),  and 
Lemma  3.4,  assuming  that  such  asymptotic  forms  exist  for  impedance  functions.  Then,  we 
demonstrate  the  existence  of  such  asymptotic  expansions  (Lemma  3.9),  by  converting  the 
Schrodinger  equation  into  an  integral  equation  (Lemma  3.5)  and  applying  the  Neumann 
series  (Lemma  3.6).  Finally,  the  statements  (A)  and  (B)  are  formulated  and  demonstrated 
by  Theorems  3.10,  3.12.  The  following  lemma  yields  the  hrst  two  terms  in  the  asymptotic 
forms  of  the  impedance  functions  p_, 


Lemma  3.4  Suppose  that  impedance  functions  p+{x,  k) ,  p-{x,k),  defined  by  (32),  (33),  are 
given  by  the  asymptotic  series 


/  /X  aAx)  aoix) 

p+{x,  k)  =  ao{x)  +  +  •  ■  •, 

p_{x,  k)  =  ho{x)  +  ^  +  ■  •  •, 

tk  [iky 


for  large  real  k.  Then, 


ao{x)  =  bo{x)  =  — -  ■  ^/lT~q{xy+Ty^^, 
p{x) 


(133) 

(134) 


(135) 


aAx)  =  —bAx)  = 


p{x)  ■  {q'{x)  +  i  fi{x))  —  2  ■  (1  -I-  q{x)  +  i  jix)  —  a^)  ■  p'{x) 


4  ■  p‘^{x)  ■  (1  -|-  q(x)  -|-  i  'y{x)  —  a^) 


(136) 


where  p,  q,  7,  a  are  defined  in  (13). 


Proof.  It  is  easily  observed  from  (82),  (83)  that  the  impedance  functions p+(a;,  —k),  p-{x,  k) 
satisfy  the  same  Riccati  differential  equation  (83).  Hence,  we  have 


afix)  =  bfix),  for  all  the  even  i, 


(137) 


afix)  =  —hi{x),  for  all  the  odd  i.  (138) 

(135),  (136),  as  well  as  other  asymptotic  coefficients  hi{x),i  =  2,3,...,  are  obtained  by 
plugging  (134)  into  (83)  and  comparing  terms  of  different  orders  of  k.  □ 


The  following  lemma  converts  the  Schrodinger  equation  (69)  into  an  integral  equation. 
Lemma  3.5  Suppose,  for  all  x  &  R  and  complex  k  y  0, 


RX 

/  +  q^r)  —  a‘^  +  iji^r)  dr, 

Jo 

(139) 

m{t,  k)  =  k), 

(140) 

^ikt 

n{t,k)=  y_{t,k), 

(141) 
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whereip-  is  defined  in  (67),  q,p,'y  E  Cq([0,  1])  sueh  that  l+q{x)  —  a‘^  >  0,  7(0;)  >  0,  p{x)  >  0. 
Suppose  further  that  T  is  a  path  defined  in  the  complex  plane  such  that 


{r  :  t  G  r,a;  G  R,t{x)  =  /  a/1  +  g(r)  —  +  i  j(t)  dr}. 


(142) 


Then,  the  differential  equation  (69)  subject  to  boundary  condition  (71)  can  be  converted  into 
integral  equations  by  the  formulae 


m  =  Fk{m)  +  (1  +  gi  +  *71  -  % 


where 


nit,  k)  =  m(t,  k) 


mm  = 


2ik 


7(r)  k)  dr, 


2ik 


ri(T){l-r‘''*''-pf(T)dT. 


(143) 

(144) 

(145) 


Proof.  Combining  (69),  (71)  with  (140),  (141),  we  observe  that  m  satishes  the  equation 

m''{t,  k)  —  2ikm'(t,  k)  =  —p{t)  m{t,  k),  (146) 

subject  to  the  initial  conditions 

m(0, /c)  =  (1  +  gi  +  *7i  -  0^)3  %  (147) 

m'(0,A;)=0.  (148) 

Multiplying  (146)  by  and  integrating  the  result  from  0  to  t,  we  have 


m'{t,  k)  =  -  p(r)  m(r,  k)  dr. 


(149) 


(143)  is  obtained  immediately  via  integrating  (149)  from  0  to  t,  and  equation  (144)  follows 
from  (149),  (143),  (144).  □ 


Observation  3.2  Since  pipr)  is  continuous  on  the  entire  complex  plane  and  zero  outside  of 
a  bounded  region,  the  functions  and  p(r)  bounded  for  all  real 

t,  T,  and  k  G  .  Therefore,  there  exists  a  real  number  ci  >  0  such  that 


\Fk\  < 


Cl_ 

\ky 


and  hence  there  exists  a  real  number  A  >  0  such  that 


(150) 


Fk\  <  1, 


(151) 


for  all  k  G  K{A),  defined  in  Section  2.1. 
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Remark  3.3  For  complex,  not  purely  real  t  and  r,  (150)  does  not  hold.  Due  to  (59)  and 
(60),  it  is  easy  to  observe  that  t,  r  are  real  if  and  only  if  the  attenuation  does  not  enter 
into  our  scattering  problem,  i.e.  7(0;)  =  0  for  all  x  &  R.  Therefore,  (150)  holds  only  in  the 
absence  of  attenuation  7.  While  our  analysis  applies  only  to  the  case  when  7  =  0,  numerical 
experiments  in  Section  5  indicate  that  our  scheme  still  works,  when  the  attenuation  is  small, 
i.e.  |7(a;)|  |1  +  q{x)  —  a^|,  for  all  x  E  R,  When  the  attenuation  is  relatively  large,  our 

scheme  does  not  work. 


Lemmas  3.6  and  3.7  yield  the  Neumann  series  for  the  integral  equation  (143)  and  an  estimate 
of  the  error  for  using  a  truncation  of  the  series. 

Lemma  3.6  Suppose  that  q,  p,^  E  Cg([0, 1]),  p>2,  q^^^\  7*^^)  are  absolutely  continuous 

for  all  X  E  R.  Suppose  further  that  T  is  a  path  defined  in  the  complex  plane  such  that 


{r  :tET,xE  R,t{x)  =  /  a/1  +  g(r)  —  a‘^  +  i  7(r)  dr}.  (152) 

Jo 

Then  for  any  integer  I  <  I  <  p,  there  exist  Oj  :  T  ^  R,  j  =  1, p  —  1,,  :  T  x  C'^  C, 

where 

d^~^aj{t) 
dti^-i 

are  bounded  and  absolutely  continuous  for  all  t  E  T,  j  =  1,  ...,p  —  1,  and 

a^{t,  k)  (154) 

is  bounded  and  absolutely  continuous  function  oft  for  all  (t,  k)  eT  x  C~^,  such  that 

R—  1  - 

mi{t,k)  =  {1  +  qi  +  i^i-a^)^  pf^  + 

i=i 

where  :  F  x  C  is  defined  by  the  formulae 


mo(t,  k)  =  0, 


(156) 


mz(t,  /c)  =  (1  +  gi  +  i7i  -  0^)4  2  Ffc(mz_i)(t,  k) 

=  (1  +  gi +i7i /  p(r)(l  -  e^*^^‘“^))mz_i(r, /c)  dr.  (157) 

Zl  K  Jq 

Proof.  We  prove  this  lemma  by  induction.  For  I  =  1,  formulae  (156),  (157)  yield 

mi(f,  A;)  =  (1  +  gi  +  *71  -  0^)2  Pi  ^  (158) 

for  all  {t,  k)  eT  X  C^,  which  is  already  in  the  form  (155)  satisfying  conditions  (153),  (154). 
For  I  >  1,  assuming  that  mi(t,k)  is  in  the  form  (155)  satisfying  conditions  (153),  (154),  we 
obtain  m/+i  using  (157): 

m«+i(t, /c)  =  (1  +  gi  +  *7i  -  0^)2  p^  2  f  p(r)(l  -  A;)  dr 

Zl  K  Jq 

=  (1  +  O'!  +  *  7i  —  0^)4  Pi  ^  +  di(t,  k)  +  Rit,  k)  +  Rit,  k)  +  R(t,  k)  (159) 
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with  Ij  :  r  X  C~^  C ,  1  <  j  <  4  defined  by  the  formulae 

h{t,  k)  =  ^  ^  riir)  +  (^)  /  aj-i{r)  dr,  (160) 

hit,  k)  =  vir)  (1  -  dr,  (161) 

Isit,  k)  =  -J2  (^)  ^  dir)  a.-i(r) 

M-i 

=  -J2hit,k),  (162) 

s=2 

liit,  k)  =  ^  ^  77(r)  a^(r)(l  -  ^t.  (163) 

Clearly,  we  only  need  to  show  that  /j  ,  1  <  j  <  4  can  be  expressed  in  the  form 

+  (164) 

j=i  V  / 

with  ttj  :  T  — ^  -R,  1  <  j  <  p  —  1  satisfying  condition  (153)  and  :  F  x  — >  C  satisfying 


condition  (154).  Obviously,  Ji  and  are  already  in  the  form  (164).  We  now  use  Lemma  2.1 
to  show  that  I2,  h  can  also  be  expanded  in  the  form  (164).  Observing  that  rjitix))  =  0  for 
all  X  ^  (0, 1),  770-2)  jg  absolutely  continuous,  and  that  1  <  j  <  /i  —  1  are  absolutely 

continuous  (due  to  the  inductive  assumption),  we  can  use  formula  (40)  to  expand  I2  and 
each  term  Jsis  =  1, ...,  p  —  1)  of  I3  as 


biit,k), 


(165) 


hit,  k)  = 


2ik 


s  pt 


^  /  1  1) 


■Z 

J=S  +  1 


2ikJ 


—rrivit)  as-lit)) 


riir)  as-iir)  e^^^d-r) 

I  \  ^^ 


2ik 


bsit,k) 


(166) 


with  bs  '■  T  X  C~^  C  uniformly  bounded  on  F  x  (see  Lemma  2.1).  Therefore,  I2  is  in 
the  form  (164)  due  to  (165),  and  h  is  of  the  form  (164)  due  to  (166),  (162).  Thus,  mi+iit,  k) 
can  indeed  be  written  in  the  form  (155)  satisfying  conditions  (153),  (154).  □ 
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Lemma  3.7  Suppose  that  the  funetions  m,  n,  :  R  x  C  are  defined  by  the 

formulae  (147),  (148),  (157)  and 


n^(t,  k)  =  m^it,  k) 


2ik 


k)  dr 


(167) 


respectively.  Then  under  the  conditions  of  the  preceding  lemma,  there  exist  positive  real 
numbers  A,  ci,  C2,  C3  such  that 


\m{t,k)  -  mf,{t,k)\  <  — 

\n{t,k)  -  n^{t,k)\  < 
for  all  {Re{t),  k)  G  R  x  K{A),  and 

I  n{t,k)  ,  ^  C3 
m{t,  k)  ~  \kfi^ 


(168) 

(169) 


(170) 


for  all  {Re{t),k)  G  [Ti,cx))  x  K{A). 

Proof.  Due  to  (150),  the  norm  of  the  integral  operator  Fj^  in  (157)  is  of  the  order  0(|/c|“^) 
for  any  k  G  C~^,  from  which  we  observe  that  there  exists  A  >  0,  such  that  (168)  is  true. 
Subtracting  (167)  from  (144),  we  obtain 


n(t,  k)  -  n^it,  k) 

=  m{t,k)  -  mf,{t,k)  -  [  pir)  k)  -  k))  dr.  (171) 

Zl  K  Jq 

Now,  the  estimate  (169)  is  a  direct  consequence  of  (171),  (168),  and  the  fact  stated  in 
Observation  3.2  that  in  (171),  the  expression 

,172) 

is  bounded  for  all  k  G  K{A),  —00  <  r  <  t  <  00  .  We  now  prove  (170)  by  showing  that  there 
exists  a  positive  number  C3  such  that 


nf,{t,k)  ^  C3 

m^{t,k)  ~  \kfi 


(173) 


for  all  {Re(t),  k)  G  [i?e(Ti),  cx))  x  K{A).  Lemma  3.6  states  that  a^(t,  k)  in  (155)  is  bounded 
and  absolutely  continuous  for  all  {t,k)  G  L  x  ,  and  aj{t)  in  (155)  is  also  independent 
of  k.  Therefore,  we  can  assume  that  the  constant  A  has  been  chosen  such  that  for  all 
{Re{t),k)  eRx  K{A), 


.7=1 
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or  equivalently, 

I'rrif.it,  k)\>^  -  {l  +  qi+i'ji-  a‘^)^  (175) 

Combining  (167)  with  (155),  we  obtain 

nf,  =  m^  +  hit,  k)  +  hit,  k)  +  hit,  k),  (176) 

with  hit,  k),  hit,  k)  dehned  by  (161),  (162),  and  hit,  k)  dehned  by  the  formula 

75(C/i^)  =  (^)  j  p(r)a^(r,A;)e2*^(*“^)dr.  (177) 

Noticing  that  p(r)  =  0  for  all  i?e(t)  >  i?e(Ti),  we  have 

hit,k)  =  biit,k),  (178) 

Jsit^k)  =  bsit,k),  (179) 

for  all  (i?e(t),  k)  G  [i?e(Ti),  cx))  x  KiA),  due  to  (165),  (166).  Consequently,  there  exists  c  >  0 
such  that 

\hit,  k)  +  hit,  k)  +  hit,  k)\  <  (180) 

for  all  (i?e(t),  k)  G  [i?e(Ti),  cx))  xKiA),  since  a^(t,  k),  bgit,  k)  are  bounded  for  all  (i?e(t),  k)  G 
[i?e(Ti),  cx)  X  KiA),  and  s  =  1, ...,  n  —  1. 

Now,  the  estimate  (173)  is  a  direct  consequence  of  (176),  (180)  and  (175).  The  esti¬ 
mate  (170)  is  a  direct  consequence  of  (173),  (168),  and  (169).  □ 

The  proof  of  the  following  lemma  is  to  that  of  Lemma  3.7,  and  is  therefore  omitted. 


Lemma  3.8  Suppose,  for  all  x  &  R  and  complex  k  ^  0, 


tix)  =  /  a/1  -|-  g(r)  —  -f  f  7(r)  dr. 


(181) 


/(«,*;)  =e-“V>+(«X),  (182) 

_ 2  /t  ^ 

git,  k)  =  ^—^h'+it,  k),  (183) 

whereh+  is  defined  in  (66),  q,p,'y  E  Co([0, 1])  such  that  l+qix)  —  a‘^  >  0,  7(0;)  >  0,  pix)  >  0. 
Then  under  the  conditions  of  the  Lemma  3. 1,  there  exist  positive  numbers  A,  d^  such  that 


9it,k)  ^  da 

fit,k)  ~  \kfi' 


(184) 


for  all  (i?e(t),  k)  G  (— cx,  0]  x  KiA). 
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Now,  we  are  ready  to  show  the  existence  of  the  asymptotic  expansion  (133),  (134) 
for  impedance  functions  p_|_,  p_  by  converting  Schrodinger  equation  into  an  integral  equa¬ 
tion  (Lemma  3.5)  and  using  the  Neumann  series  (Lemma  3.6). 

Lemma  3.9  Suppose  that  impedance  functions  p+{x,k)  p-{x,k)  are  defined  by  (32),  (33), 
for  all  X  E  R,  k  E  C.  Then,  there  exist  asymptotic  series  expansions  for  p+  and  p_;  that 
is,  there  exist  sequences  of  complex  functions  a  =  {oj  :  R  ^  C},  and  b  =  {bi  :  R  ^  C}, 
i  =  0, 1,2, ...  such  that  and  p-  are  asymptotic  given  by  the  series 

/  /X  aiix)  aoix) 

p+{x,k)  =  ao{x)  +  —-^  +  —^  +  ---  (185) 

tk  (iky 

p.{x,k)=ho{x)  +  ^-^+^-^  +  --;  (186) 

as  \k\  oo. 

Proof.  Combining  (74)  with  (168),  (169),  (140),  (141),  (155),  and  (176),  we  obtain  (186). 
(185)  is  derived  similarly.  □ 

Theorems  3.10,  3.12  concern  the  statements  (A)  and  (B)  outlined  in  the  beginning  of 
Section  3. 

Theorem  3.10  Suppose  that  g,p,  7  G  Co([0, 1]),  1  -|-  q{x)  —  >  0,  7(0;)  >  0,  p{x)  >  0  for 

all  X  E  R  and  q" ,  p" ,  7"  are  absolutely  continuous.  Suppose  further  that 

D  =  {(x,  k)\x  E  R,  Im{k)  >  0}.  (187) 


Then 

(a)  0+  and  0.  are  continuous  functions  of  {x,  k)  and  analytic  functions  of  k  for  allx  E  R 
and  k  E  C ; 

(b)  p+  and  p_  are  continuous  functions  of  {x,  k)  and  analytic  functions  of  k  in  D; 

(c)  For  all  {x,  k)  E  D, 


p+{x,  k)  =  ——  ■  ^Jl  +  q{x)  +*7 
p{x) 


1  p{x)  ■  {q\x)  -I-  i  '^'{x))  —  2  ■  (1  -I-  q{x)  +  i  7(0;)  —  a^)  ■  p'{x) 


i  k 


4  ■  p‘^{x)  ■  (1  -|-  q{x)  +  i  'j{x)  —  ay 


+  0{k-y,  (188) 


P-{x,  k)  =  ——  ■  A/r+”g(xyTT7^^^ 


p{x 


1 

i  k 


p{x)  ■  {q'{x)  +  i  7'(a;))  -  2  ■  (1  -h  q{x)  +  i  'y(x) 


ay  ■  p'(x) 


4  ■  p^(x)  ■  (1  -|-  q(x)  -f  i  ^(x)  —  a^) 


+  0(k-y,  (189) 
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Proof.  We  only  give  the  proof  for  p_  since  the  proof  for  (j)^,  is  very  similar.  We 
introduce  two  auxiliary  functions  01,02^  R  x  C~^  ^  C  hy  the  formulae 

^{x,  k)  =  (p-{x,  k)  —  1,  (190) 

(f){x,k)  =  —  _q,2  (XgX) 

p{x)  Pi 

and  combining  (190),  (191)  with  (13)  and  initial  conditions  (34),  (35),  we  obtain  the  linear 
hrst  order  ODEs 


0'(a;,  k)  =  pix)  ^(x,  k)  —  i  k  +  qi  +  i  yi  —  a‘^ - , 

Pi 


(p'{x,  k)  = 


p{x] 


subject  to  the  initial  conditions 


-(1  +  q{x)  +i'y{x)  -  Q;^)(0(a;,  k)  +  1), 

0(0, /c)  =  0, 

0(0,  k)  =  0. 


(192) 

(193) 

(194) 

(195) 


According  to  Lemma  2.4,  0,  0  are  continuous  functions  of  x  and  analytic  functions  of  k 
for  all  a;  G  i?  and  k  &  C,  from  which  part  (a)  follows  immediately.  Similarly,  we  obtain 
part  (b)  by  combining  part  (a)  with  (33)  and  the  fact  that  (f)-{x,k)  ^  0  for  all  {x,k)  G  D 
(see  Lemma  3.3).  The  expansion  (189)  follows  immediately  from  Lemmas  3.9  and  3.4.  □ 

Corollary  3.11  Under  the  conditions  of  the  preceding  theorem,  there  exist  positive  number 
Cl,  C2  such  that 

pk  ff  p+{T,k)  p(t)  drl 


<  Cl, 

(196) 

<  C2, 

(197) 

I  g*  If  V-  U,k)  p{t)  (It 

for  all  t,x  E  [0, 1],  k  E  R,  or  for  all  0  <  t  <  x  <  1,  k  E  C~^ . 

Proof.  Due  to  parts  (b)  and  (c)  of  Theorem  3.10,  the  real  part  of  the  functions 

(198) 

(199) 

where  C3  and  C4  does  not  depend  on  t  or  a;  for  t,x  E  [0, 1],  k  E  R,  or  for  all  0  <  t  <  a;  <  1, 
from  which  (196)  (197)  follow  immediately.  □ 


Re  k  I  p+{T,k)dTj  <  C3, 
Re(ik  f  p-{T,k)d'j\  <04, 
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Theorem  3.12  Suppose  that  q,''^,p  G  c™([0, 1]),  m  >  2,  7^”^)  are  absolutely 

continuous  and  1  +  q{x)  —  >  0,  7(0;)  >  0,  p{x)  >  0  for  all  x  E  R.  Then  there  exists  a 

positive  real  number  a  such  that 

\p+{x,-k) -p_{x,k)\  <  (200) 

for  all  {x,  k)  E  R  X  C~^ 

Proof.  Due  to  (74),  (75),  we  obtain 

pUx.-i)  -P-U.fc)  =  Vl+pW+«7W-a^-3^  ■  (201) 

Combining  Lemmas  3.7,  (34),  3.8  and  (35)  yields  that  (201)  is  true  for  all  x  ^  (0, 1).  In  order 
to  prove  the  theorem  for  x  E  (0, 1),  we  observe  that  p+{x,  —k)  and  P-{x,  k)  obey  the  same 
Riccati  equation  (83)  due  to  (82)  and  (83).  The  difference  s{x,k)  =  p+{x,—k)  —p^{x,k) 
satisfies  the  ODE 

s'{x,  k)  =  ik  p{x)  {p+{x,  —k)  +p_(a;,  k))  s{x,  k). 

Clearly,  the  solution  to  (202)  is 

s{x,k)  =  e*^/o"(p+C-D+p-(i,D)pWdts(o,  A;). 

(196)  and  (197)  show  that  there  exists  constant  &  >  0  such  that 

\^^ik  {p+{t-k)+p-{t,k))  pit)  dt\^  <  1) 

for  all  {x,  k)  E  [0, 1]  x  R.  Due  to  (85),  (74),  and  Lemma  3.8,  there  exists  a  positive  number 
c  such  that  for  all  k  E  R, 

|s(0,/i;)|  =  \p+{0,-k)  -p-{0,k)\  =  \p+{0,-k)  -  ^  ^  _c^  ^^05) 

Pi  \k\ 

Now,  (200)  for  x  E  (0, 1)  follows  immediately  from  (203),  (205).  Thus,  we  have  (200)  for  all 
{x,k)ERxC+.  □ 


(202) 

(203) 

(204) 


3.3  Trace  Formula 

In  this  section,  we  prove  Theorem  3.13,  which  is  the  principal  analytical  tool  of  this  paper. 
Theorem  3.13  describes  what  are  known  as  the  trace  formulae  for  the  impedance  functions 
P+,  p-  in  the  context  of  varying  density,  speed  of  propagation,  and  attenuation. 

Theorem  3.13  (Trace  Formula).  Suppose  that  q,p,'y  E  c™([0, 1]),  m  >  2,  q^^\ 
pim)  absolutely  continuous  and  1  +  q{x)  —  >  0,  7(0;)  >  0,  p{x)  >  0  for  all  x  E  R. 

Then, 

p{x)  ■  {q'{x)  +  i  'y'{x))  —  2  ■  p\x)  ■  (1  +  q{x)  +  i  j(x)  —  a^) 

2 

=  —  {I  +  q{x)  +  i'y{x)  —  a'^)  p‘^{x)  /  {p+{x,  k)  —  p-{x,  k))  dk.  (206) 

J  —  00 
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Moreover,  there  exists  a  positive  number  c  such  that 
\p{x)  ■  {q{x)  +  i  y(x))  -  2  ■  p'(x)  ■  (1  +  q(x)  +  i  j(x)  -  a^) 

- -(l  +  q(x) +ij(x) -a^)p^(x)  [  {p+{x,k)  -  p_{x,k))  dk\  <  J  (207) 

for  all  X  E  R,  a  >  0. 

Proof.  Due  to  part  (C)  of  Theorem  3.10,  there  exists  c  >  0  such  that 


\{p+{x,  k)  -p-{x,  k))- 

1  p(x)-(q'(x)+iy(x))-2-(l  +  q(x)+ij(x)-a^)-p'(x)  c 

ik  2  ■  p'^{x)  ■  [1  +  q{x)  +  i'^{x)  —  a^)  ~  1^1^ 

for  all  {x,  k)  E  R  X  C~^.  Denoting  by  T  the  upper  half  circle  of  radius  A,  with  clockwise 
orientation,  in  the  complex  /c-plane,  i.e.. 


T  =  {k\kEC+,\k\  =A}, 

and  noting  that  —  p_  is  an  analytical  function  of  k  E  C~^,  we  obtain 


{p+{x,  k)  —  p-(x,  k))  dk  =  /  (P+  {x,  k)  —  p-{x,  k))  dk. 


>-A 


t'X 


Substituting  (208)  into  (210),  we  have 


—  {1  +  q{x)  +  i'-){x)  —  a^)  p^{x)  ■  /  {p+{x,k)  —  p-{x,k))  dk  = 


TT 


'-A 


p{x)  ■  {q'{x)  +  i  ^\x))  —  2  ■  p'{x)  ■  (1  +  q{x)  +  i  j(x)  —  a^)  +  0{k~^) 
from  which  (206)  follows  immediately.  In  order  to  prove  (207),  we  rewrite  (206)  as 
p{x)  ■  {q{x)  +  i  'y'{x))  —  2  ■  p{x)  ■  (1  +  q{x)  +  i  'y{x)  —  a^)  = 


71 


{1  +  q{x)  +  i'y{x)  —  a  )  p  {x)  ■  /  {p+{x,  k)  —  p_{x,  k))  dk  +  I{a) 


with  I  (a)  given  by  the  formula 
2 

I{a)  =  —  (1  +  q{x)  Riji^x)  —  a^)  p‘^{x)  ■ 

71 


—a  /*cxD> 

+ 

— cxD  Ja 


Due  to  the  symmetry  of  the  integrals  in  (213),  we  have 


/(a)  =  -  (1  +  q{x)  +i'^{x)  -  a  )  p  (x) 

TT 


—a  roQ"^ 
+ 

—  oo  o  a 


{p+{x,  k)  —  p_{x,k))  dk. 


{p+{x,  —k)  —  p-{x,  k))  dk, 


and  using  (200),  we  obtain  a  constant  c  such  that 


|/(a)|  < 


c 


a 


(m— 1)  ’ 


from  which  (207)  follows  immediately. 


(209) 

(210) 

(211) 

(212) 

(213) 

(214) 

(215) 

□ 
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Remark  3.4  As  there  are  two  unknowns  in  (206),  (view  q{x)  +  i'y{x)  as  one  complex  un¬ 
known,  p{x)  as  the  other),  at  least  two  a’s  are  needed  in  order  to  solve  p',  q' ,  r'  in  (206).  In 
particular,  we  obtain  the  trace  formulae 


p\x) 


Re(F(Q;i))  —  Re(F(Q;2)) 
2  («?  -  ai) 


(216) 


q'{x) 

with 


.  w  2 - ^  (Re(F(ai))  (1  -h  q{x)  -  aj)  -  Re(F(a2))  (1  +  q{x)  -  a^)) 

P[X)  tCli  Oi2) 

_  Im(F(ai))  ■  (aj  -  a|)  7(0;)  ■  (Re(F(ai))  -  Re(F(a2))) 

^  “  Pix)  (a?  -  ai) 

F(a)  =  - 

IT 


{1  +  q{x)  +  i'-^{x)  —  a^)  p‘^{x)  /  {p+{x,  k)  —  p^{x,  k))  dk, 


(217) 

(218) 

(219) 


for  the  case  of  two  a’s.  Multiple  (more  than  two)  choices  of  a’s  would  lead  to  an  overdeter¬ 
mined  complex  linear  system,  and  thus  can  be  used  to  control  the  effects  of  noise. 


4  The  Algorithm 

This  section  describes  the  algorithm  of  the  present  paper  and  gives  details  about  its  imple¬ 
mentation  and  computational  costs. 

4.1  Description  of  the  Algorithm 

In  this  section,  we  describe  a  reconstruction  algorithm  for  the  scalar  Helmholtz  equation  in 
layered  acoustic  media 

4>±{x,  k)  -  ^  ■  (/)±(x,  k)  +  k^  ■  (1  +  q(x)  +  i  ■  j(x)  -  a^)  ■  0±(a;,  k)  =  0,  (220) 

p{x) 

subject  to  the  initial  conditions 

(j)+{x,  k)  =  72-^2 X  ^  ^221) 

0_(a;,  k)  =  e-*^Vi+'?i+*7i-«2^  for  all  a;  <  0.  (222) 

In  (220)  -  (222),  a;  is  a  real  number,  /c  is  a  complex  number  in  the  upper  half  plane,  a  is  the 
sine  of  the  angle  of  incidence  with  respect  to  the  normal  to  the  interface  of  layers,  0+  and 
0_  are  the  scalar  fields  associated  with  right-going  and  left-going  waves,  respectively;  the 
parameters  to  be  recovered  in  this  algorithm  are  the  density  p,  potential  g,  and  attennation 
7  of  the  layered  media.  We  assnme  p,  g,  7  G  c™([0, 1]),  i.e.,  p,  g,  7  have  m  continnous 
derivatives  everywhere,  and  are  dehned  by  eqnations  (6)  -  (11). 

As  discussed  in  Sections  2  and  3,  in  order  to  reconstruct  parameters  p,  g,  7,  we  consider 
a  system  of  integro-differential  eqnations 

4(x.  k)  =  t  p(n  .  k)  -  (223) 
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p'_  (x,  k)  =  ik  p{x)  ■  (  pt  {x,  k)  — 
p\x)  = 


1  +  q{x)  +  i  ^{x)  —  a" 


p^{x) 

Re(F(Q;i))  —  Re(F(Q;2)) 


2(«f 


Q^o 


q'{x)  = 


(Re(F(Q;i))  (1  +  q{x)  -  a^)  -  Re(F(a2))  (1  +  q(x)  -  a^)) 


p(x)  (af  -  al) 

lm(F(ai))  ■  (af  -  a\)  +  7(0;)  ■  (Re(F(Q;i))  -  Re(F(Q;2))) 


7  [X)  = 


p{x)  (af  —  al 


with 


F{a)  =  —  {1  +  q{x)  +  iji^x)  —  p‘^{x)  /  {p+{x,  k)  —  p-{x,  k))  dk, 


TT 


p+{0,k)  =po{k) 
a/1  +  (?i  +  *71 


subject  to  the  initial  conditions 

p-i0,k)  = 

Pi 

P(0)  =  Pi, 

g(0)  =  qi, 

7(0)  =  7i. 

In  (223),  (224),  the  impedance  functions  p+,p_  :  {R,C~^) 

(j)'^{x,  k) 


P+{.x,  k)  =  - 
P-{x,  k)  = 


i  k  p{x)  0+(a;,  k)  ’ 
(j)'_{x,  k) 


(224) 

(225) 

(226) 

(227) 

(228) 

(229) 

(230) 

(231) 

(232) 

(233) 

C  are  dehned  by  the  formulae 

(234) 


—i  k  p{x)  (x,  k) '  (235) 

(223)  and  (224)  are  Riccati  equations  obtained  directly  from  the  Helmholtz  equation  (220) 
and  the  definitions  of  impedance  functions  (234),  (235);  equations  (225),  (226),  (227)  are 
known  as  trace  formulae,  connecting  the  Fourier  components  of  the  solutions  of  the  Helmholtz 
equation  to  the  parameters  of  the  scattering  objects  to  be  recovered. 

Our  reconstruction  algorithm  for  the  inverse  scattering  problem  in  layered  acoustic  media 
amounts  to  solving  numerically  a  self-contained  set  of  ODEs,  i.e.,  (223)  -  (227),  subject  to 
the  initial  conditions  (229)  -  (233).  In  this  paper,  the  ODE  solver  from  [8]  is  used. 

As  we  shall  see  in  Section  5,  for  sufficiently  large  a,  the  system  of  ODEs  (223)  -  (227)  has  a 
unique  solution  for  all  a;  G  [0, 1],  and  this  solution  is  stable  with  respect  to  small  perturbations 
of  the  initial  data  po{k).  The  inversion  algorithm  is  (m  —  l)*^-order  convergent  for  all  three 
parameters  p,  g,  7  to  be  recovered,  where  m  is  the  smoothness  of  p,  g,  and  7. 


4.2  Implementation 

In  implementing  the  algorithm  stated  above,  the  integral 


{p+{x,  k)  —  p-{x,  k))  dk 


(236) 


in  (228)  is  approximated  by  the  trapezoidal  rnle  T„,  i.e., 


Tn{p+{x,k) -p-{x,k))  =  h  {p+{x,kj) -p-{x,kj)) 


j=-M+l 


+  -P-{x,-a))  +  {p+{x,a)  -p_{x,-))), 


(237) 


with  h  =  a/M,  kj  =  j  h,  j  =  —M, M.  Thus,  the  system  of  integro-differential  equations 
(223)  -  (227)  subject  to  initial  conditions  (229)  -  (233)  is  converted  into  a  system  of  8M  +  7 
ODEs 

//  •/  t  \  (  2t  l+q{,x) +i-i{x)  - 

P+{x,  kj)  =  -i  kj  p{x)  ■  I  p+(a;,  kj) - -  1  (238) 


p_{x,  kj)  =  i  kj  p{x)  ■  ^t{x,  kj)  — 


1  +  q{x)  +  i  'y{x)  — 


p'{x)  = 


Re(F(Q;i))  —  Re(F(Q;2)) 
2  (a?  -  al) 


(239) 


(240) 


q'{x)  = 


p{x)  (aj  -  al) 


(Re(F(Q;i))  (1  +  q{x)  -  aj)  -  Re(F(Q;2))  (1  +  qix)  -  a^))  (241) 


f  _  Im(F(ai))  ■  {af  -  a^)  +  7(x)  ■  (Re(F(ai))  -  Re(F(a2))) 
^  “  p(x)  (a?  -  ai) 


(242) 


F(a)  =  -  (1  +  g(x)  +ij(x)  -  a^)p^(x)  ■  T2m+i(p+(x,  k)  -p_(x,k)),  (243) 

71 


and  subject  to  the  initial  conditions 


p+{0,kj)  =po{kj 


p-{0,kj)  = 


_  a/1  +  gi  +i7i  - 

Pi 

P(0)  =  Pi, 

g(0)  =  gi, 

7(0)  =  7i. 


(244) 

(245) 

(246) 

(247) 

(248) 
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Remark  4.1  In  all  the  numerical  examples  in  Section  5,  the  values  of  the  initial  impedance 
functions  Po{kj),  j  =  — M,  required  for  the  reconstruction  scheme,  are  provided  by 

solving  forward  scattering  problems,  namely,  AM  +  1  independent  ODEs 


0"(a;,  kj)  -  ■  0'(a;,  kj)  +  kj  ■  {1  +  q{x)  +  i  ■  'y{x)  -  ■  (j){x,  kj)  =  0, 

p[x) 

subject  to  the  boundary  conditions 


(249) 


0+(a;,  kj)  =  for  all  a;  >  1,  (250) 

for  kj  =  j  ■  jj,  j  =  — M, ...,  M  and  a  =  ai,a2.  Again,  we  used  the  ODE  solver  in  [8]. 
Remark  4.2  Due  to  Observation  2.3,  for  all  x,k  E  R, 

(251) 

(252) 


p+{x,k)  =p+{x,-k), 


P-{x,  k)  =  p-{x,  —k), 

thus,  the  integral  /“^(p+(a;,  k)  —  p-{x,  k))  dk  in  (228)  is  equal  to 


2-  /  'Re{pj^{x,k)  —  p-{x,k))  dk. 


(253) 


Therefore,  the  dimensions  of  the  system  of  ODEs  we  consider  (see  equations  (238)  -  (242) 
is  reduced  to  AM  +  7  from  8M  +  7. 


Remark  4.3  According  to  Lemma  3.2  and  Remark  3.1,  impedance  functions  p^{x,k), 
p^{x,k)  are  continuous  in  the  vicinity  of  A;  =  0.  This  allowed  us  to  use  Lagrange  inter¬ 
polation  to  get  the  values  of  p+{x,  k),  p-{x,  k)  ed.  k  =  0. 


4.3  Complexity  Analysis 

The  time  cost  of  the  inverse  scheme  is  of  the  order  0{Nk  ■  N^),  where  Nj^  is  the  number  of 
measurements  in  the  frequency  domain,  and  is  the  number  of  nodes  in  the  space  domain, 
since  the  computational  cost  for  the  ODE  solver  we  use  is  proportional  to  the  dimension 
of  the  ODE  system  [Nk  in  our  case)  and  the  number  of  discretization  points  in  the  space 
domain  {Nz). 

Further,  the  storage  requirements  of  the  algorithm  are  also  determined  by  and  N^, 
and  is  of  the  form 

S  =  0{K-Nk)+0{Nz),  (254) 

where  A'  is  a  constant  determined  by  the  precision  required  by  the  ODE  solver  in  [8].  For 
single  precision,  K  =  22;  for  double  precision,  K  =  60. 
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5  Numerical  Examples 

The  algorithm  of  Section  4  has  been  implemented  in  Fortran  77  in  double  precision.  In  this 
section,  we  illustrate  the  performance  of  the  scheme  as  applied  to  several  different  classes  of 
scattering  objects,  from  Gaussian  to  discontinuous  staircase-shaped  ones.  The  experiments 
were  carried  out  on  a  2.8GHz  Pentium  D  desktop  with  2  Gb  of  RAM  and  an  L2  cache 
of  1  Mb.  The  calculations  reported  in  Examples  1  and  2  were  carried  out  with  a  requested 
accuracy  of  10“^®  in  the  ODE  solver;  the  calculations  reported  in  Examples  3-7  were  carried 
out  with  a  requested  accuracy  of  10“^. 

In  Examples  1-3,  the  scatterers  satisfy  the  smoothness  conditions  of  Theorem  3.3.  In 
Examples  4  -  4.4,  the  scatterers  violate  the  smoothness  conditions  mildly,  as  the  scatterers 
are  continuous  but  their  derivatives  are  not  continuous.  In  Examples  5-6,  the  scatterers 
strongly  violate  the  smoothness  conditions,  as  those  scatterers  are  discontinuous.  In  Example 
7,  we  performed  a  crude  test  of  stability  of  the  algorithm  by  truncating  the  input  scattering 
data  p+(— 1,  /c)  and  using  it  in  the  reconstruction  algorithm.  The  headings  of  the  Tables  are 
defined  as  follows: 

a  is  the  largest  frequency  used  in  the  algorithm; 

hk  is  the  step  size  in  the  discretization  of  frequency; 

Nx  is  the  number  of  discretization  points  in  [—1, 1]; 

Ep,  Eg,  E"^  are  the  relative  norms  of  error  of  p,  q,  7; 

E^,  E^,  E^  are  the  relative  maximum  norms  of  error  of  p,  q,  7; 

tcpu  is  the  GPU  time  required  in  seconds. 

Example  1  ;  In  this  example,  we  reconstruct  scattering  parameters  p,  q,  and  7  of  the 
Gaussian  distribution  given  by  the  formulae 

p{x)  =  1000  +  500  ■  (255) 

g(a;)  =  (256) 

7(0;)  =  0.01  +  0.01-6-^°^'.  (257) 

This  is  an  example  of  scatterer  whose  p,  q,  and  7  are  in  in  the  interval  [—1,1]  np  to 
double  precision.  Tables  1  and  2  illustrate  the  numerical  behavior  of  the  reconstruction 
algorithm,  and  Figure  1  contains  graphs  of  the  exact  and  the  recovered  p,  g,  7  G  (they 
are  almost  indistinguishable  in  the  graph)  and  the  input  impedance  function  p^{—l,k).  In 
this  example,  the  algorithm  converges  extremely  rapidly  as  we  would  expect. 

Example  2  :  In  this  example,  we  reconstruct  a  more  complicated  scattering  object  given 
by  the  formulae 


p{x) 

=  1000  +  1000  ■  6“^°^'  ■  cos(30a;) 

^bx  _ 

^bx  _|_  ^  —  bx  ’ 

(258) 

^5  X 

^—5  X 

q{x)=e  ^°"  ■cos(20a;)- 

(259) 

^bx  y-,—bx 

7(0;) 

=  0.01  +  0.004  ■  6"^°^'  ■  sin(10a;) 

6  —  6 
gSa;  _j_  g— 5x  ’ 

(260) 
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(1.1)  Reconstruction  of  p 


Figure  1: 


(1.2)  Reconstruction  of  1  +  g 
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(1.4)  Real  Part  of  p+(— 1,  A;) 


of  Example  1  with  a  =  50 


Figure  2:  Reconstruction  of  Example  2  with  a  =  100 


This  is  another  example  of  a  smooth  scatterer  in  the  interval  of  [-1,1],  i.e.,  p,  q,  7,  G 
C'“[— 1,1].  Tables  3  and  4  illustrate  the  numerical  behavior  of  the  algorithm;  graphs  con¬ 
taining  the  exact  and  the  reconstructed  p,  q,  7  as  well  as  the  input  impedance  functions 
p+(— 1,/c)  are  included  in  Figure  2.  The  algorithm  also  yields  superalgebraic  convergence, 
although  not  as  fast  as  that  of  Example  1  due  to  the  highly  oscillating  character  of  the 
scatterer  and  the  impedance  functions,  as  can  be  seen  in  Figure  2. 
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Table  1:  CPU  time  and  accuracies  for  Example  1 


Table  2:  CPU  time  and  accuracies  for  Example  1 


Table  3:  CPU  time  and  accuracies  for  Example  2 


Table  4:  CPU  time  and  accuracies  for  Example  2 

a  hk  Nx _ _ _ E^ _ tcpu 

50  0.05  1000  2.30E-02  1.03E-01  1.22E-02  5.9E+01 

100  0.1  1000  1.62E-03  1.05E-02  8.04E-03  6.0E+01 

100  0.05  1000  9.99E-05  4.14E-04  1.57E-04  1.4E+02 

100  0.025  2000  1.06E-04  4.55E-04  1.71E-04  5.0E+02 

200  0.1  2000  1.04E-03  6.83E-03  5.23E-03  2.3E+02 

200  0.05  2000  9.83E-08  5.99E-07  3.55E-07  5.0E+02 

200  0.05  4000  9.28E-08  5.62E-07  3.26E-07  9.0E+02 
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Example  3  :  In  this  example,  we  reconstruct  a  scatterer  defined  by  the  formulae 


t  =  {x  +  1)  ■  Tl 


(261) 


,  22  6 
p{x)  =  1000  +  100  ■  (1  -  cos(4t))  -  —  (1  -  cos(5t))  +  —  (1  -  cos(7t))  , 

25  49 

1215  7 

q{x)  =  0.4  ■  (  (1  -  cos(3t))  -  (1  -  cos(llt))  +  ^  (1  “  cos(12t))  )  , 

16  5 

7(0;)  =  0.01  +  0.003  ■  (  (1  —  cos(2t)) - (1  —  cos(3t))  H - (1  —  cos(4t))  )  . 

21  28 


(262) 

(263) 

(264) 


The  scatterer  is  a  Cg-function  in  R  with  support  in  the  interval  [—1,1].  The  performance  of 
the  algorithm  is  demonstrated  in  Tables  5,  6  and  Figures  3,  4.  As  we  can  see  from  those 
tables,  the  convergence  of  the  algorithm  is  actually  better  than  our  prediction  of  4*^-order 
convergence. 
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0.016 

0.014 

0.012 

0.01 

0.008 

0.006 

0.004 


-1  -0.5  0  0.5  1 

(3.3)  Reconstruction  of  7 


Figure  3:  Reconstruction  of  Example  3  with  a  =  50 
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Table  5:  CPU  time  and  accuracies  for  Example  3 


a 

hk 

N. 

tcpu 

25 

0.025 

MHiin 

7.95E-02 

2.58E-01 

6.78E-02 

1.7E+01 

50 

0.025 

1.86E-02 

6.57E-02 

1.84E-02 

4.1E+01 

100 

0.025 

1.87E-04 

4.72E-04 

3.34E-04 

3.5E+02 

200 

0.0125 

8.50E-07 

2.55E-06 

9.45E-07 

2.7E+03 

Table  6:  CPU  time  and 

.  accuracies 

For  Example  3 

a 

hk 

N, 

ET 

ET 

E^ 

tcpu 

25 

0.025 

MHIIII 

2.44E-01 

1.24E+00 

4.40E-01 

1.7E+01 

50 

0.025 

5.24E-02 

2.20E-01 

1.05E-01 

4.1E+01 

100 

0.025 

5.18E-04 

1.56E-03 

1.12E-03 

3.5E+02 

200 

0.0125 

2.12E-06 

7.46E-06 

3.86E-06 

2.7E+03 
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Example  4  :  In  this  example,  we  construct  a  scatterer  with  discontinuous  derivatives 
supported  on  [—1,1],  dehned  by  the  formulae 


p{x)  =  1000  +  500  ■  sin(7a;). 

(265) 

Zx  _  ^ 

g(a;)  =  0.2-cos(30a;)- 

(266) 

^\lx  _  ^—\lx 

-i{x)  =  0.01  +  0.004  ■  cos(20a;)  ■  ^ 

(267) 

Tables  7  and  8  illustrate  the  numerical  behavior  of  the  algorithm;  Figure  5  demonstrate 
the  exact  and  the  reconstructed  p,  q,  7,  and  the  input  impedance  function  p_|_(— 1, /c).  The 
algorithm  is  not  convergent  in  this  case,  although  the  input  impedance  function  p+(— 1, /c) 
behaves  properly.  Further  investigations  (see  Figure  6)  show  that,  as  we  move  with  the  ODE 
solver  towards  the  right  boundary  of  the  scattering  structure,  both  the  impedance  function 
p+{x,  k)  and  the  integrand  p+(a;,  k)  —p-{x,  k)  in  the  trace  formulae  (206)  become  extremely 
oscillatory  and  blow  up  as  k  increases.  This  phenomena  is  closely  related  to  the  interaction 
between  non-smooth  behavior  of  p,  q  and  the  accumulated  effect  of  attenuation.  Example 
4.1  and  Example  4.2  explore  this  phenomenon  in  more  detail.. 

Example  4.1  ;  This  example  uses  the  same  p  and  q  as  in  Example  4,  but  with  zero 
attenuation;  thus  we  have 


p{x)  =  1000  +  500  ■  sin(7a;). 

(268) 

^Zx  g— 3a; 

g(a;)  =  0.2-cos(30a;)-^3^_^^_3^. 

(269) 

7(0;)  =  0. 

(270) 

Table  9  illustrates  the  numerical  behavior  of  the  reconstruction  algorithm,  and  Figure  7  con¬ 
tains  graphs  of  the  exact  and  the  recovered  p,  q.  The  algorithm  exhibits  linear  convergence. 
Unlike  Example  4,  the  integrand  p+{x,  k)  —p-{x,  k)  in  the  trace  formulae  (206),  and  thus  in 
the  ODE  system,  is  not  increasing  with  k  as  the  wave  travels  through  the  scattering  object 
(see  Figure  8). 

Example  4.2  :  This  example  uses  the  same  p  and  q  as  in  Example  4  and  the  same  7  as 
in  Example  3  (cg-function);  thus  we  have 


t  =  (x  -f  1)  ■  TT 


p{x)  =  1000  -I-  500  ■  sin(7a;). 


q{x)  =  0.2  ■  cos(30a;) 


g3 X  I  g  Zx 


j{x)  =  0.01  +  0.003- 


16  5 

(1  -  cos(2t))  -  ^  (1  -  cos(3t))  +  ^  (1  -  cos(4t)) 


(271) 

(272) 

(273) 

(274) 


Tables  10,  11  illustrate  the  numerical  behavior  of  the  reconstruction  algorithm. 

Example  4.3:  This  example  uses  the  same  q  and  7  as  in  Example  4,  but  with  a  constant 
p;  thus  we  have 

p{x)  =  1000,  (275) 
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Table  7:  CPU  time  and 

.  accuracies 

for  Example  4 

a 

hk 

N. 

K 

E'i 

tcpu 

25 

0.025 

MHiin 

1.74E-02 

4.10E-02 

1.73E-01 

1.7E+01 

50 

0.025 

3.22E-02 

1.88E-01 

1.66E-01 

4.0E+01 

100 

0.025 

2.11E-02 

1.30E-01 

8.29E-02 

1.7E+02 

200 

0.025 

4000 

1.71E-02 

1.13E-01 

1.62E-01 

6.7E+02 

Table  8:  CPU  time  and  accuracies 

for  Example  4 

a 

hk 

N. 

ET 

ET 

E^ 

tcpu 

25 

0.025 

MHiin 

4.97E-02 

1.04E-01 

4.89E-01 

1.7E+01 

50 

0.025 

■nM 

1.27E-01 

7.17E-01 

6.67E-01 

4.0E+01 

100 

0.025 

2000 

7.74E-02 

4.58E-01 

3.40E-01 

1.7E+02 

200 

0.025 

4000 

8.94E-02 

5.44E-01 

7.49E-01 

6.7E+02 

Table  9:  CPU  time  and  accuracies  for  Example  4.1 


a 

hk 

N. 

ET 

ET 

tcpu 

25 

|HI^ 

WHIIII 

1.92E-02 

6.32E-02 

5.43E-02 

1.81E-01 

8.2E+00 

50 

■nM 

3.15E-02 

1.85E-01 

1.22E-01 

6.91E-01 

1.7E+01 

100 

0.025 

2000 

1.54E-02 

9.59E-02 

5.87E-02 

3.20E-01 

8.0E+01 

200 

0.05 

2000 

5.98E-03 

3.68E-02 

2.35E-02 

1.17E-01 

8.0E+01 

400 

0.05 

4000 

3.33E-03 

2.06E-02 

1.31E-02 

6.90E-02 

3.5E+02 

800 

0.05 

4000 

1.97E-03 

1.22E-02 

7.70E-03 

4.15E-02 

7.0E+02 

Table 

-0:  CPU  time  and 

accuracies 

for  Example  4.2 

a 

hk 

N. 

E'^k 

tcpu 

25 

0.025 

1.74E-02 

3.98E-02 

1.19E-02 

1.7E+01 

50 

0.025 

2000 

3.26E-02 

1.90E-01 

3.68E-02 

8.1E+01 

100 

0.025 

2000 

2.34E-02 

1.43E-01 

4.01E-02 

1.7E+02 

200 

0.05 

4000 

3.86E-02 

2.44E-01 

1.83E-01 

3.5E+02 

Table  11:  CPU  time  and 

accuracies 

for  Example  4.2 

a 

hk 

N. 

ET 

ET 

ET 

tcpu 

25 

4.93E-02 

1.03E-01 

5.33E-02 

1.7E+01 

50 

2000 

1.28E-01 

7.25E-01 

l.OlE-01 

8.1E+01 

100 

0.025 

2000 

8.46E-02 

5.14E-01 

9.85E-02 

1.7E+02 

200 

0.05 

4000 

2.50E-01 

1.63E+00 

9.21E-01 

3.5E+02 

40 


0.0004 


0.0003 


0.0002 
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0  20  40  60  80  100120140160180200  20  40  60  80  100  120  140  160  180  200 

(6.1)  Re(p+(-l, /c))  (6.2)  Re  (p+(-l, /c)  -  p_(-l, /c)) 


0.0022 

0.002 

0.0018 

0.0016 

0.0014 

0.0012 


0.0008 

0.0006 

0.0004 

0.0002 


n  nnnoc 


0  . . .  -0.00015  . . . 

0  20  40  60  80  100120140160180200  '  20  40  60  80  100  120  140  160  180  200 

(6.3)  Re(p+(0,/c))  (6.4)  Re  (p+(0, /c)  -  p_(0, /c)) 


0.0009 

0.00088 

0.00086 

0.00084 

0.00082 

0.0008 

0.00078 

0.00076 

0.00074 

0.00072 

0.0007 

0.00068 


0  20  40  60  80  100120140160180200 

(6.5)  Re(p+(1,  A;)) 


0.00015 

1e-04 

5e-05 


-5e-05 

-0.0001 

-0.00015 

-0.0002 

-0.00025 


20  40  60  80  100  120  140  160  180  200 

(6.6)  Re(p+(1,  A;)  -p_{l,k)) 


Figure  6:  Reconstruction  of  Example  4  with  a  =  200 


Figure  7:  Reconstruction  of  Example  4.1  with  a  =  200 


0  tX/ 

g(a;)  =  0.2-cos(30a;)-^3^^^_3^. 

(276) 

pllx  p  —  llx 

7(x)  =  0,01  +0.004.cos(20x).^„^^^_„.. 

(277) 

Table  12  illustrates  the  numerical  behavior  of  the  reconstruction  algorithm.  Figure  9  contains 
the  exact  and  the  reconstructed  q  and  7,  and  the  input  impedance  function  p+(—l,  k).  Linear 
convergence  is  observed  for  both  q  and  7. 

Example  4.4:  This  example  uses  the  same  discontinuous  7  as  in  Example  4,  but  with 
constant  p,  g;  thus  we  have 

p{x)  =  1000,  (278) 

q{x)  =  0,  (279) 

7(a;)  =  0.01  +  0.004  ■  cos(20  x)  ■  (280) 

Table  13  illustrates  the  linear  convergence  of  the  reconstruction  algorithm,  as  we  can  expect 
given  Example  4.3.  Figure  10  contains  the  exact  and  the  reconstructed  7,  and  the  input 
impedance  function  p_|_(— 1,  /c). 
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0.009 

0.008 

0.007 

0.006 

0.005 

0.004 

0.003 

0.002 

0.001 

0 


0  20  40  60  80  100120140160180200 

(8.1)  Re(p+(-l, /c)) 


0.0006 

0.0005 

0.0004 

0.0003 

0.0002 

1e-04 

0 

-0.0001 

-0.0002 

-0.0003 

-0.0004 


20  40  60  80  100  120  140  160  180  200 

(8.2)  Re{p+{-l,k) -p_{-l,k)) 


0.003 

0.0025 

0.002 

0.0015 

0.001 

0.0005 

0 


0  20  40  60  80  100120140160180200 

(8.3)  Re(p+(0,/c)) 


0.0002 

0.00015 

1e-04 

5e-05 

0 

-5e-05 

-0.0001 

-0.00015 

-0.0002 


20  40  60  80  100  120  140  160  180  200 

(8.4)  Re(p+(0,/c)  -p_(0,/c)) 


0.000771 

0.00077 

0.000769 

0.000768 

0.000767 

0.000766 

0.000765 

0.000764 

0.000763 

0.000762 

0.000761 


0  20  40  60  80100120140160180200 

(8.5)  Re(p+(1,  A;)) 


0.00025 

0.0002 

0.00015 

1e-04 

5e-05 

0 

-5e-05 

-0.0001 

-0.00015 

-0.0002 


20  40  60  80  100  120  140  160  180  200 

(8.6)  Re{p+{l,k)-p_{l,k)) 


Figure  8:  Reconstruction  of  Example  4.1  with  a  =  200 
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Table  12:  CPU  time  and  a 


a 

hk 

N, 

E': 

^7 

25 

0.05 

MHIIII 

7.73E-02 

1.35E- 

50 

0.1 

4.46E-02 

1.20E- 

100 

0.1 

2.22E-02 

7.67E- 

200 

0.1 

4000 

1.15E-02 

1.74E- 

(9.1)  Reconstruction  of  1  +  g 


0  20  40  60  80 

(9.3)  Real  Par 


Figure  9:  Reconstruction 


ncuracies  for  Example  4.3 


ET 

tcpu 

01 

2.06E-01 

4.67E-01 

8.1E+00 

01 

1.24E-01 

4.63E-01 

8.5E+00 

02 

6.61E-02 

3.41E-01 

3.3E+01 

02 

3.53E-02 

6.12E-02 

1.6E+02 

0.015 

0.014 

0.013 

0.012 

0.011 

0.01 

0.009 

0.008 

0.007 

0.006 

0.005 


-1  -0.5  0  0.5 

(9.2)  Reconstruction  of  7 


T 


_ I _ I _ I _ I _ I _ I 

100120140160180200 

t  of  p+(— 1,  k) 


Example  4.3  with  a  =  200 


Table  13:  CPU  time  and  accuracies  for  Example  4.4 


a 

hk 

N. 

tcpu 

25 

0.05 

4.71E-02 

1.08E-01 

8.1E+00 

50 

0.1 

Iniiiil 

2.26E-01 

4.51E-02 

8.0E+00 

100 

0.2 

2000 

l.llE-02 

2.06E-02 

8.2E+01 

200 

0.2 

8000 

5.74E-03 

1.02E-02 

6.7E+02 

0.001025 

0.00102 

0.001015 

0.00101 

0.001005 

0.001 

0.000995 

0.00099 

0.000985 

0.00098 

0.000975 


(10.1)  Reconstruction  of  7 


0  20  40  60  80100120140160180200 

(10.2)  Real  Part  of  p+(— 1,  k) 


Figure  10:  Reconstruction  of  Example  4.4  with  a  =  200 
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Example  5  :  In  this  example,  we  reconstruct  a  scatterer  defined  by  the  formulae 


1 

r  1000 

X  G 

(— CX),  —0.5] 

p(a;)  =  \ 

1200 

X  G 

T 

0 

0 

(281) 

I 

[  1000 

X  G 

(0.5,  cx) 

r  0 

X  G  ( 

1 

8 

1 

0 

q{x)  = 

{ 0.4 

X  G  ( 

-0.5,  0.5] 

(282) 

1 0 

X  G  (0.5,  cx)) 

f 

0.010 

X  G 

(— X,  —0.5] 

7(a;)  =  1 

0.015 

X  G 

T 

0 

0 

(283) 

1 

0.010 

X  G 

(0.5,  x) 

In  this  example,  p,  q,  7  are  discontinuous.  Tables  14,  15  illustrates  the  numerical  behavior  of 
the  reconstruction  algorithm,  and  Figure  11  contains  graphs  of  the  exact  and  the  recovered 
p,  q,  7  and  the  input  impedance  function  p^{—l,k).  The  algorithm  does  not  converge  as 
we  can  expect  from  Example  4;  The  integrand  p+{x,  k)  —  P-{x,  k)  is  highly  oscillatory  and 
blows  up  as  k  increase  for  all  a;  G  i?  (see  Figure  12). 

Example  5.1  In  this  example,  we  reconstruct  a  scatterer  whose  attenuation  7  is  discon¬ 
tinuous  and  density  p,  potential  q  are  constant,  given  by  the  formulae 


p{x)  =  1000, 
q{x) 


(284) 


l{x)  = 


=  0,  (285) 

0.010  a;  e  (-CX), -0.5] 

0.015  a;  e  (-0.5, 0.5]  .  (286) 

0.010  a;  e  (0.5,  00) 

The  numerical  results  are  demonstrated  in  Table  16;  the  exact  and  the  reconstructed  7,  as 
well  as  the  input  impedance  function  p+(— 1,  k)  are  contained  in  Figure  13.  The  convergence 
of  the  algorithm  is  of  the  order  O(^),  where  a  is  the  largest  frequency  in  the  reconstruction 
algorithm.  The  convergence  of  the  algorithm  for  discontinuous  7  is  quite  similar  to  that  of 
discontinuous  p,  as  reported  in  [3]. 
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Table  14:  CPU  time  and  accuracies  for  Example  5 


a 

hk 

N. 

E'i 

tcpu 

25 

0.05 

1.18E-02 

2.58E-02 

3.08E-02 

4.2E+00 

50 

0.05 

7.94E-03 

2.02E-02 

2.27E-02 

1.7E+01 

100 

0.025 

4000 

1.36E-02 

2.92E-02 

2.78E-02 

1.7E+02 

200 

0.025 

1.28E-02 

2.76E-02 

2.53E-02 

1.3E+03 

Table  15:  CPU  time  and  accuracies  for  Example  5 


a 

hk 

N. 

ET 

ET 

tcpu 

25 

0.05 

7.75E-02 

1.91E-01 

2.39E-01 

4.2E+00 

50 

0.05 

6.43E-02 

1.78E-01 

2.13E-01 

1.7E+01 

100 

0.025 

4000 

1.25E-01 

3.52E-01 

4.93E-01 

1.7E+02 

200 

0.025 

Dll 

1.14E-01 

3.44E-01 

4.77E-01 

1.3E+03 

Table  16:  CPU  time  and  accuracies  for  Example  5.1 


a 

hk 

N. 

E'^ 

tcpu 

25 

0.05 

■nniil 

3.13E-02 

2.36E-01 

8.3E+00 

50 

0.1 

2.20E-02 

2.37E-01 

8.7E+00 

100 

0.2 

1.56E-02 

2.35E-01 

8.5E+00 

200 

0.2 

1.27E-02 

2.26E-01 

6.9E+01 
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(11.1)  Reconstruction  of  p 


Figure  11: 


(11.2)  Reconstruction  of  1  +  g 


0.00105 

0.00104 

0.00103 

0.00102 

0.00101 

0.001 

0.00099 

0.00098 

0.00097 

0.00096 

0.00095 


1  0  20  40  60  80  100120140160180200 

(11.4)  Real  Part  of  p+(— 1,  k) 


of  Example  5  with  a  =  50 


(13.1)  Reconstruction  of  7  (13.2)  Re  (p+(— 1,  k)) 

Figure  13:  Reconstruction  of  Example  5.1  with  a  =  200 
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Example  6  :  Here,  we  reconstruct  a  staircase-shaped  scatterer  defined  by  the  formulae 


p{x) 


q{x) 


l{.x) 


1050 

X  G 

(— cx),  —0.8 

1150 

X  G 

T 

o 

bo 

1 

o 

1250 

X  G 

(-0.4,  0.0] 

1350 

X  G 

(0.0,  0.4] 

1300 

X  G 

(0.4,  0.8] 

1200 

X  G 

(0.8,  cxd) 

'  0 

X  G 

(— CX),  —0.8] 

0.1 

X  G 

(-0.8, -0.6] 

0.2 

X  G 

(-0.6, -0.2] 

0.3 

X  G 

(-0.2,  0.2] 

0.2 

X  G 

;o.2,o.8] 

0 

\ 

X  G 

(0.8,  cx) 

0.01 

X  G 

(—X,  —0.8] 

0.012 

X  G 

T 

o 

bo 

1 

o 

0.01 

X  G 

(-0.6, -0.2] 

0.008 

X  G 

(-0.2,  0.2] 

0.007 

X  G 

(0.2,  0.6] 

0.008 

X  G 

(0.6,  0.8] 

0.009 

X  G 

(0.8,  x) 

(287) 


(288) 


(289) 


The  example  is  similar  to  the  preceding  one,  but  the  shape  of  the  scatterer  is  more  compli¬ 
cated.  The  numerical  results  are  shown  in  Table  17,  18  and  Figures  14,  15,  16,  and  17.  The 
algorithm  does  not  converge  in  this  situation. 
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-0.5  0  0.5 

(14.1)  Reconstruction  of  p 


Figure  14: 


(14.4)  Real  Part  of  p+(— 1,  k) 


of  Example  6  with  a  =  25 


-1  -0.5  0  0.5 

(15.1)  Reconstruction  of  p 


Figure  15: 


(15.4)  Real  Part  of  p+(— 1,  k) 


of  Example  6  with  a  =  100 


(16.1)  Reconstruction  of  p 


I  I  I 


(16.3)  Reconstruction  of  7 


Figure  16: 


(16.2)  Reconstruction  of  1  +  g 


(16.4)  Real  Part  of  p+(— 1,  k) 


of  Example  6  with  a  =  200 


0.0012 

0.00115 

0.0011 

0.00105 

0.001 

0.00095 

0.0009 

0.00085 

0.0008 

0.00075 


(17.1)  Re  (p+(-l,  k))  (17.2)  Re  (p+(-l,  k)  -  p_(-l,  k)) 


(17.3)  Re  {p+{0,k)) 


(17.4)  Re(p+(0,/c)  -p_(0,/c)) 


(17.5)  Re(p+(1, /c)) 


(17.6)  Re{p+{l,k)  -p.{l,k)) 


Figure  17:  Reconstruction  of  Example  6  with  a  =  200 
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Table  17:  CPU  time  and  accuracies  for  Example  6 


a 

hk 

N, 

tcpu 

25 

0.05 

■ifiM 

1.81E-02 

3.77E-02 

2.83E-02 

4.2E+00 

50 

0.05 

1000 

2.55E-02 

5.50E-02 

2.25E-02 

1.8E+01 

100 

0.1 

1000 

2.10E-02 

4.53E-02 

2.39E-02 

1.8E+01 

200 

0.05 

4000 

3.77E-02 

7.85E-02 

5.06E-02 

3.5E+02 

Table  18: 


CP 


U  time  and  accuracies  for  Example  6 


a 

hk 

N. 

ET 

ET 

ET 

tcpu 

25 

0.05 

4.99E-02 

1.05E-01 

1.79E-01 

4.2E+00 

50 

0.05 

1000 

5.69E-02 

1.23E-01 

1.49E-01 

1.8E+01 

100 

0.1 

1000 

7.22E-02 

1.78E-01 

1.95E-01 

1.8E+01 

200 

0.05 

4000 

2.72E-01 

6.57E-01 

2.71E-01 

3.5E+02 
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(18.1)  Reconstruction  of  p 


Figure  18:  Reconstruction  of  Example  1  with  D=3 


Example  7.  In  this  example,  we  investigate  the  sensitivity  of  the  reconstruction  to  per¬ 
turbations  of  the  initial  data.  We  perturb  the  initial  data  for  the  algorithm  by  truncating  it 
after  a  specified  number  D  of  decimal  digits  (both  the  real  and  the  imaginary  parts).  Thus, 
the  maximum  relative  error  is  of  the  order  10^“^,  e.g.,  when  the  number  1.999  is  truncated 
after  D  =  2  digits,  the  result  is  1.9. 

Tables  19,  20  demonstrates  the  numerical  results  of  the  reconstruction  of  Example  1  with 
various  truncations  of  the  input  data.  In  each  case,  a,  hk,  is  chosen  properly  so  that 
the  error  caused  by  the  inaccuracy  of  the  input  data  is  much  larger  than  that  caused  by  the 
insufficient  a,  hk  or  iV^,.  Figures  18,  19,  20,  21  demonstrate  the  exact  and  the  reconstructed 
p,  q  and  7  for  different  orders  of  perturbations. 
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Table  19:  Accuracies  for  Example  1  with  Truncated  Data 


D 

a 

hk 

N. 

3 

50 

0.0125 

4000 

1.52E-01 

3.72E-01 

3.65E-01 

3 

100 

0.0125 

8000 

1.47E-01 

3.57E-01 

3.50E-01 

4 

50 

0.025 

2000 

4.97E-02 

9.52E-02 

9.47E-02 

4 

100 

0.025 

4000 

4.99E-02 

9.53E-02 

9.48E-02 

5 

50 

0.025 

2000 

l.OOE-02 

2.07E-02 

2.08E-02 

5 

100 

0.025 

4000 

l.OlE-02 

2.07E-02 

2.09E-02 

6 

50 

0.00625 

8000 

5.92E-04 

1.20E-03 

1.19E-03 

6 

100 

0.00625 

5.94E-04 

1.20E-03 

1.20E-03 

Table  20:  Accuracies  for  Example  1  with  Truncated  Data 


D 

a 

hk 

N. 

ET 

ET 

E^ 

3 

50 

0.0125 

8.22E-01 

2.32E+00 

2.31E+00 

3 

100 

0.0125 

5.09E-01 

1.27E+00 

1.27E+00 

4 

50 

0.025 

2000 

1.16E-01 

2.29E+00 

2.28E+00 

4 

100 

0.025 

4000 

1.37E-01 

2.56E+00 

2.57E+00 

5 

50 

0.025 

2000 

2.02E-02 

3.97E-02 

4.00E-02 

5 

100 

0.025 

4000 

2.09E-02 

4.23E-02 

4.23E-02 

6 

50 

0.00625 

8000 

1.56E-03 

2.91E-03 

2.90E-03 

6 

100 

0.00625 

16000 

1.60E-03 

3.01E-03 

2.99E-03 
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0.022 
0.02 
0.018 
0.016 
0.014 
0.012 
0.01 
0.008 
0.006 

-1  -0.5  0  0.5  1 

(19.3)  Reconstruction  of  7 

Figure  19:  Reconstruction  of  Example  1  with  I 


60 


The  following  observations  can  be  made  from  the  tables  above,  and  from  a  wide  range  of 
numerical  tests  performed  by  us. 

1.  For  scatterers  satisfying  the  conditions  of  Theorem  3.3  (Example  1,  2,  3),  the  numerical 
algorithm  of  Section  4  displays  convergence  of  order  m  —  1,  where  m  is  the  smoothness  of 
the  scatterer;  the  CPU  time  required  is  proportional  to  N^.  ■  N^,  where  and  are  the 
numbers  of  discretization  points  in  frequency  and  space  domain,  respectively. 

2.  For  general  scatterers  violating  the  conditions  of  Theorem  3.3  mildly,  the  algorithm 
does  not  converge.  However,  the  algorithm  exhibits  linear  convergence  (see  Example  4.1, 
4.3,  4.4)  for  the  following  two  particular  categories  of  scatterers  violating  the  conditions  of 
Theorem  3.3  mildly, 

A.  p,  q  are  continuous  but  their  derivatives  are  not,  and  7  is  absent 

B.  q,  7  are  continuous  but  their  derivatives  are  not,  and  p  is  a  constant. 

3.  When  the  scatterer  is  discontinuous  (Example  5,  5.1,  6),  the  algorithm  produces 
results  demonstrated  in  Figures  11  -  17.  The  oscillatory  behavior  near  the  discontinuities 
is  the  well-known  Gibbs  phenomenon.  In  general,  the  algorithm  is  not  convergent  for  such 
scatterers.  However,  for  scatterers  of  the  following  categories 

A.  q  is  discontinuous,  p  is  a  constant,  7  is  absent, 

B.  7  is  discontinuous,  p  and  q  are  constants, 

the  convergence  of  the  algorithm  is  of  the  order  O(^),  where  a  is  the  largest  frequency. 

4.  When  the  initial  data  is  perturbed  (Example  7),  the  error  of  the  reconstruction  is 
proportional  to  the  magnitude  of  the  perturbation,  and  the  proportion  coefficient  is  1. 


6  Conclusions 

In  this  paper,  we  construct  numerical  algorithms  for  the  solution  of  inverse  scattering  prob¬ 
lems  in  layered  acoustic  media  in  three  dimensions.  The  speed  c  of  propagation  of  sound,  the 
density  p,  and  the  attenuation  7  are  the  three  parameters  reconstructed  by  the  algorithm. 
The  computational  complexity  of  the  algorithm  is  0{Nk  ■  N^),  where  W,  are  the  number 
of  discretization  points  in  the  frequency  and  space  domains,  respectively. 

The  inverse  scattering  schemes  we  construct  can  be  viewed  as  an  extension  of  [3] ,  and  are 
based  on  a  collection  of  trace  formulae,  which  connect  the  Fourier  components  of  the  solutions 
of  the  Helmholtz  equation  to  the  parameters  being  recovered.  Under  the  assumption  of  mild 
attenuation  and  a  smoothly  varying  medium,  our  inverse  scattering  algorithms  require  only 
a  few  measurements,  in  the  sense  that  given  a  medium  whose  parameters  c,  p,  and  7  have 
m  >  1  continuous  derivatives,  and  data  measured  for  all  frequencies  u  in  the  interval  [—a,  a], 
the  error  of  the  reconstruction  decays  as  as  a  ^  00.  In  this  respect,  our  algorithm 

is  similar  to  the  Fourier  transform. 

Acknowledgments.  The  authors  would  like  to  thank  M.  Tygert  for  constructive  dis¬ 
cussions. 
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